How to find my ode45 equation in specific h

Question

esat gulhan 2020-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582449-how-to-find-my-ode45-equation-in-specific-h

编辑： Alan Stevens 2020-8-21

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-8-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582449-how-to-find-my-ode45-equation-in-specific-h#answer_482993

编辑：Alan Stevens 2020-8-21

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Alan Stevens 2020-8-21

在 MATLAB Online 中打开

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 2020-8-21

编辑：Alan Stevens 2020-8-21

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;

请先登录，再进行评论。