How to solve this equation?

Hey guys

How can I solve this simple ode equation with matlab??

[Wan Ji]

Hi 金成勋 my friend,

This is a first-order nonlinear ode with power exponent other than 1, far more complicated than you ever imagine.

Answer 1

[Star Strider]

One approach —

syms h(t) t Y 
Dh = diff(h);
Eqn = 25*pi*Dh == pi*(4*0.0254)^2 * sqrt(2*9.8*h*2*(101325/876)+Dh^2)+90/3600

Eqn(t) = 

Eqn = isolate(Eqn,Dh)

Eqn = 

[VF,Subs] = odeToVectorField(Eqn)

VF = 

Subs = 

VF = simplify(VF, 500)

VF = 

hfcn = matlabFunction(VF, 'Vars',{t,Y})

hfcn = function_handle with value:
    @(t,Y)[(pi.*8.112963841460668e+32+sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)]

Then use the appropriate numerical ODE solver (most likely ode15s) to integrate it.

.

[Wan Ji]

formidable, ami!

[성훈 김]

i'm really sorry but how can i use appropriate ODE solver??

I'm a really beginner in matlab.

[Star Strider]

No worries!

Use ‘hfcn’ as it exists, with the appropriate time vector and initial conditions:

hfcn = @(t,Y)[(pi.*8.112963841460668e+32+sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)];

tspan = [0 100];
ic = 0;
[t,y] = ode15s(hfcn, tspan, ic);

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

With the sqrt call, this is likely half of the actual solution, so a negative solution likely also exists. I leave that to you to determine (and to code).

Make appropriate changes to get the result you want. (I ran this on my computer in order to get the ‘hfcn’ expression, in order to avoid running the earllier code again here.)

.

[성훈 김]

i'm sorry to bother you..ㅠㅠ

I want to find the expression h(t-2) with a negative slope.

Can you help me please?

[Star Strider]

The equation itself does not have any parameters that can be estimated that would give a negative slope.

The sqrt term has both positive and negative roots, so change the sign of that term to get the negative square roots:

hfcnp = @(t,Y)[(pi.*8.112963841460668e+32+sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)];
hfcnn = @(t,Y)[(pi.*8.112963841460668e+32-sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)];
%                                       ↑ ← HERE

tspan = [0 10];
ic = 0;
[tp,yp] = ode15s(hfcnp, tspan, ic);
[tn,yn] = ode15s(hfcnn, tspan, ic);

figure
yyaxis left
plot(tp, yp)
ylabel('h(t) +Root')
yyaxis right
plot(tn, yn)
grid
xlabel('t')
ylabel('h(t) -Root')
legend('Positive Root','Negative Root', 'Location','SE')

I doubt that it has a negative slope anywhere.

The only way to force that would be to negate the derivative:

% hfcnp = @(t,Y)-[(pi.*8.112963841460668e+32+sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)];
% hfcnn = @(t,Y)-[(pi.*8.112963841460668e+32-sqrt(7.3e+1).*sqrt(Y(1).*-1.129628243491514e+37+pi.^2.*Y(1).*6.713376166760685e+42+1.480615901066572e+33).*3.201062735323997e+12)./(pi.^2.*8.112963841460668e+35-1.365130281111817e+30)];
% %                                       ↑ ← HERE
% 
% tspan = [0 10];
% ic = 0;
% [tp,yp] = ode15s(hfcnp, tspan, ic);
% [tn,yn] = ode15s(hfcnn, tspan, ic);
% 
% figure
% yyaxis left
% plot(tp, yp)
% ylabel('h(t) +Root')
% yyaxis right
% plot(tn, yn)
% grid
% xlabel('t')
% ylabel('h(t) -Root')
% legend('Positive Root','Negative Root', 'Location','SE')

I cannot run that (the commented-out code) here because it times out, and even takes an extraordinarily long time on my computer when I run it offlline, so it may not have s solution.

.