Can't numerically solve the integral

4 visualizaciones (últimos 30 días)
Wesley Rivers
Wesley Rivers el 14 de Nov. de 2019
Comentada: Walter Roberson el 15 de Nov. de 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;

Respuesta aceptada

Walter Roberson
Walter Roberson el 14 de Nov. de 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 comentarios
Wesley Rivers
Wesley Rivers el 14 de Nov. de 2019
Thank you again. Now I know to watch for if my function is returning a zero in the denominator. Do you have any tips for altering my bounds so I don't divide by zero but still return a similar answer?
Walter Roberson
Walter Roberson el 15 de Nov. de 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

Iniciar sesión para comentar.

Más respuestas (0)

Community Treasure Hunt

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

Start Hunting!

Translated by