Can't numerically solve the integral
7 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Walter Roberson
2019-11-14
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 个评论
Walter Roberson
2019-11-15
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!