Can't numerically solve the integral
8 views (last 30 days)
Show older comments
Wesley Rivers
on 14 Nov 2019
Commented: Walter Roberson
on 15 Nov 2019
I have the following code and I can't get the code to solve for J. It outputs an answer but no matter what I try it won't actuallly solve my integral. I have tried using vpa, and double. I've tried changing int to integral, and using the @(x) symbol. I need the variables J and denomJ to output a number. Any help would be great. Thank you in advance.
b = 25;
y(x) = 15*((0.2969*sqrt(x/25))-(0.1260*(x/25))-(0.3516*(x/25)^2)+(0.2843*(x/25)^3)-(0.1015*(x/25)^4)); %needs to be in mm
A = 2 * int(y(x),x,0,b); %area in mm
Am = A/(1000^2); %area in m
t = 0.5;
dy_dx = diff(y(x),x); %derivative of y(x) with respect to x
v = sqrt(1+(dy_dx)^2);
denomJ = (2/t)*int(v,0,b); %denomiator of J
J = (4*Am^2)/denomJ;
0 Comments
Accepted Answer
Walter Roberson
on 14 Nov 2019
syms x
b = 25;
y(x) = 15*((0.2969*sqrt(x/25))-(0.1260*(x/25))-(0.3516*(x/25)^2)+(0.2843*(x/25)^3)-(0.1015*(x/25)^4)); %needs to be in mm
A = 2 * int(y,x,0,b); %area in mm
Am = A/(1000^2); %area in m
t = 0.5;
dy_dx = diff(y,x); %derivative of y(x) with respect to x
v = sqrt(1+(dy_dx)^2);
denomJ = (2/t)*vpaintegral(v,0,b); %denomiator of J
J = (4*Am^2)/denomJ;
You cannot do it by normal means because one of the terms for v includes a division by sqrt(x/25) which is infinite at the lower bound of 0.
4 Comments
Walter Roberson
on 15 Nov 2019
Using a lower bound of eps() would probably be good enough for your purposes, but you could go as low as 1e-37 without having double(int(v,1e-37,b)) bomb out.
>> double(vpaintegral(v,0,1e-37,'reltol',1e-20))
ans =
2.79523154926173e-19
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!