Can't numerically solve the integral

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;

Answer 1

[Walter Roberson]

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.

[Wesley Rivers]

Thank you. I don't know why I didn't think that having a lower bound of 0 would effect my code. For future reference how did putting "vpa intergral" affect my code so that it would output a number. I guess what I'm trying to understand is how did you know where to put the vpa so that it would work. I tried putting vpa in other places but it didn't work.

[Walter Roberson]

how did you know where to put the vpa

Experience. I saw that the int() call was returning unevaluated, and that double() was refusing to evaluate it, so I tried vpaintegral() and it worked. It was after that that I looked and saw the division by 0.

[Wesley Rivers]

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]

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