Trying to solve a definite integral but wont give numerical solution

1. Q=150;
2. k=1.49;
3. n=0.02;
4. So=0.0048;
5. b=3;
6. g=32.2;
7. x1=0;
8. y1=3; %x2=200;
9. h=20;
10. x2=200;
11. syms y;
12. f_xy=(So-((Q*n)/(k*b*(y^(5/3))))^2)/(1-((Q^2)/((b^2)*g*y^3)));
13. func=int(f_xy,y,[x1 x2]);
14. vfunc=vpa(func);
15. y=y1+vfunc;
16. ys=string(y);
17. yf=str2double(ys);

this is what shows up in the varaibles for ys= "int((10000/(22201*y^(10/3)) - 3/625)/(12500/(161*y^3) - 1), y, 0, 200) + 3"

Also the above equation is not how it should look with variables plugged in.

Answer 1

[Torsten]

Q=150;
k=1.49;
n=0.02;
So=0.0048;
b=3;
g=32.2;
x1=0; 
y1=3; %x2=200;
h=20;
x2=200;
f_xy=@(y)(So-((Q*n)./(k*b*(y.^(5/3)))).^2)./(1-((Q^2)./((b^2)*g*y.^3)));
vfunc=integral(f_xy,x1,x2);
y=y1+vfunc

[Heather Holzer]

Thank you very much!!!

[Heather Holzer]

Im getting this warning though:

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound

on error is 1.4e-02. The integral may not exist, or it may be difficult to

approximate numerically to the requested accuracy.

> In integralCalc/iterateScalarValued (line 372)

In integralCalc/vadapt (line 132)

In integralCalc (line 75)

In integral (line 88)

In open_channel (line 13)

[Torsten]

There seem to be two singular points within the interval of integration, one at y=0 and one at y where 1-((Q^2)./((b^2)*g*y.^3)) =0. You should check whether your integral exists at all.

[Heather Holzer]

ok thank you. What if i wanted to do the same integral from x1 to x2 but x2=[20,20,200] and get a table of y values per x-value? I need to do this so that i can plot these values.

[Heather Holzer]

is this close?

sym y

x2=[20,20,200];

f_xy=(So-((Q*n)./(k*b*(y.^(5/3)))).^2)./(1-((Q^2)./((b^2)*g*y.^3)));

vfunc=int(f_xy,y,[x1 x2]);

y=y1+vfunc;

ys=table(y);

[Torsten]

Q=150;
k=1.49;
n=0.02;
So=0.0048;
b=3;
g=32.2;
x1=0; 
y1=3; %x2=200;
h=20;
x2=[20,40,60];
f_xy=@(y)(So-((Q*n)./(k*b*(y.^(5/3)))).^2)./(1-((Q^2)./((b^2)*g*y.^3)));
vfunc = zeros(numel(x2,1));
for i=1:numel(x2)
  vfunc(i) = integral(f_xy,x1,x2(i));
end
y = y1+vfunc

But - as already said - you should check whether the integral exists at all.

[Heather Holzer]

I will. Thank you so so very much for your excellence!