How to find my ode45 equation in specific h

syms D g H Do
tspan = [0 120];
mgiren=0
Do=3;
D=2/10;
h0=h; 
g=9.81;
y0 = 2;
ySol= ode45(@(t,h)(mgiren-(pi()*D^2/4*(2*g*h)^0.5)/(pi()*Do^2/4)), tspan, y0)
for t=linspace(0,100,11)
fprintf('%15.5f',t,deval(ySol,t)),;fprintf('\n')
end

My differantial code is here, dt/dh=(mgiren-(pi()*D^2/4*(2*g*h)^0.5)/(pi()*Do^2/4)),h(0)=2, h(tx)=1, how can i find tx, is there anyway to find tx, i can find it with interpolate but is there any easier way.

Answer 1

[Alan Stevens]

Something like this perhaps:

tspan = 0:120;
h0=2; 
[t, h] = ode45(@rate, tspan, h0);
tx = interp1(h,t,1);
plot(t,h,tx,1,'o'), grid
text(tx+5,1,['tx = ' num2str(tx)])
xlabel('t'),ylabel('h')

function dhdt = rate(~, h)

        mgiren=0;
        Do=3;
        D=2/10;
        g=9.81;
        dhdt = (mgiren-(pi*D^2/4*(2*g*h)^0.5)/(pi*Do^2/4));

end

Ok, I guess this still uses interpolation!

You could use fzero to find it, but, for this curve I think interpolation is far simpler.

[esat gulhan]

Actually two ways are good for this. but there is a problem when i try to reach 0. interp(h,t,0) gives error, same the fzero gives error too.

[Alan Stevens]

Tht's because you can't get to zero with the data specified. There is an analytical solution to your equation, which is most easily expressed with t as a function of h - see below. You'll notice there is a logarithmic term, which tries to calculate log(0) when both mgiren and h are zero.

tspan = 0:120;
h0=2; 

mgiren=0;
Do=3;
D=2/10;
g=9.81;
a = mgiren/(pi*Do^2/4);
b = (D/Do)^2*sqrt(2*g);

T = @(h) -2*(b*sqrt(h) + a*log(a-b*sqrt(h)))/b^2 ...
        + 2*(b*sqrt(h0) + a*log(a-b*sqrt(h0)))/b^2;
    
h = h0:-0.01:0.01;
t = T(h);

plot(t,h), grid
xlabel('t'), ylabel('h')

[Alan Stevens]

Hmm. Thinking further, the log term is zeroed all the way through because it's multiplied by a (which is zero)., so, basically, you just have a square root relationship (when a is zero).

You get NaN if you try T(0).

T(h = 0) can be found from 2*sqrt(h0)/b;