Variable in function as well as integral boundary

1 view (last 30 days)
Hi,
I have an integral I want to compute and plot for different angles theta. 0-90 degrees.
For every angle I want to compute the integral of the function from 0 up to the angle theta.
theta = 1:90;
A = 0.4*10^-19;
R = 4215 * theta.^(-1.35)*10^-9;
D_0 = 0.3;
D_c = R.*(1-cosd(theta));
D = (D_0 + D_c)*10^-9;
b = 0.2*10^-6;
Fun = @(theta) (A ./ (6 .* pi .* (D_0 + R .* ( 1-cosd(theta))).^3)) .* b .* R;
F_vdW1 = integral(Fun, 0, theta);
Error using integral
Limits of integration must be double or single scalars.
Can someone please help with what I am doing wrong?
Thanks

Accepted Answer

Torsten
Torsten on 22 Oct 2022
Edited: Torsten on 22 Oct 2022
Is it this what you want ?
If not, please clarify how the theta in the integrand and in the upper bound of the integral should be interpreted. Mathematically, your notation to use the same name for both is wrong.
theta = 1:90;
A = 0.4*10^-19;
R = 4215 * theta.^(-1.35)*10^-9;
D_0 = 0.3;
D_c = R.*(1-cosd(theta));
D = (D_0 + D_c)*10^-9;
b = 0.2*10^-6;
Fun = (A ./ (6 .* pi .* (D_0 + R .* ( 1-cosd(theta))).^3)) .* b .* R;
F_vdW1 = cumtrapz(theta,Fun);
F_vdW1 = 1×90
1.0e-30 * 0 0.0461 0.0666 0.0793 0.0881 0.0948 0.1002 0.1046 0.1083 0.1115 0.1143 0.1167 0.1189 0.1209 0.1227 0.1243 0.1258 0.1272 0.1285 0.1297 0.1308 0.1319 0.1329 0.1338 0.1347 0.1355 0.1363 0.1371 0.1378 0.1385
plot(theta,F_vdW1)
  3 Comments
Torsten
Torsten on 22 Oct 2022
Edited: Torsten on 22 Oct 2022
Are you sure that "phi" should enter the equation for R in degrees and not in radians ?
And according to your graphics, only the 6 and not the complete expression 6*pi*[D0+R*(1-cos(phi))]^3 would appear in the denominator of the integrand.
A = 0.4*10^-19;
R = @(phi) 4215 * phi.^(-1.35)*10^-9;
D_0 = 0.3;
D_c = @(phi) R(phi).*(1-cosd(phi));
D = @(phi) (D_0 + D_c(phi))*10^-9;
b = 0.2*10^-6;
Fun = @(phi)(A ./ (6 .* pi .* (D_0 + R(phi) .* ( 1-cosd(phi))).^3)) .* b .* R(phi);
theta = 1:360;
F_vdW1 = arrayfun(@(theta)integral(Fun,0,theta),theta)
F_vdW1 = 1×360
1.0e-28 * 0.6801 0.5336 0.4630 0.4187 0.3872 0.3633 0.3442 0.3285 0.3152 0.3038 0.2938 0.2850 0.2771 0.2700 0.2636 0.2577 0.2523 0.2473 0.2427 0.2384 0.2343 0.2305 0.2270 0.2236 0.2204 0.2174 0.2146 0.2119 0.2093 0.2068
plot(theta,F_vdW1)

Sign in to comment.

More Answers (1)

Bruno Luong
Bruno Luong on 22 Oct 2022
Edited: Bruno Luong on 22 Oct 2022
I have no idea if the code correspondons to the formula; I just modify your code to make it work on array
theta = 1:90;
A = 0.4*10^-19;
R = @(theta) 4215 * theta.^(-1.35)*10^-9;
D_0 = 0.3;
D_c = @(theta) R.*(1-cosd(theta));
b = 0.2*10^-6;
Fun = @(theta) (A ./ (6 .* pi .* (D_0 + R(theta) .* ( 1-cosd(theta))).^3)) .* b .* R(theta);
F_vdW1 = arrayfun(@(onetheta) integral(Fun, 0, onetheta), theta)
F_vdW1 = 1×90
1.0e-28 * 0.6801 0.5336 0.4630 0.4187 0.3872 0.3633 0.3442 0.3285 0.3152 0.3038 0.2938 0.2850 0.2771 0.2700 0.2636 0.2577 0.2523 0.2473 0.2427 0.2384 0.2343 0.2305 0.2270 0.2236 0.2204 0.2174 0.2146 0.2119 0.2093 0.2068

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!