Vtgyjfy

I am trying to numerically solve the following ordinary differential equation that I encountered in one article:
$$ y''' left( x right) =-0.5cdot y left( x right) cdot y'' left( x right) -0.05 $$

This equation has the following boundary conditions:

• At $ x=0$: $y left( x right) =y' left( x right) =0 $
  


• At $ x=infty$: $y' left( x right) =1$, $y'' left( x right) =0 $
  

In the article it is written that we can rewrite our equation in the following way
begin{align}
y' left( x right) &=f left( x right) \
f' left( x right) &=g left( x right) \
g' left( x right) &=-0.5cdot y left( x right) cdot g left( x right) -0.05
end{align}

with the following boundary conditions

• At $x=0$: $y left( x right) =f left( x right) =0 $
  


• At $x=infty$: $f left( x right) =1$, $g left( x right) =0 $
  

and that we need to assume g(x=0) as a first guess and solve our equation with the Runge Kutta fourth order method (https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods).

The problem is that I do not know how to do that in my case since 4-th order Runge-Kutta method is normally used to solve first-order ODE.

An equilibrium of the demanded type seems improbable due to the constant $-0.05$. If indeed for $xapproxinfty$ the solution were to show $y''(x)approx 0$, one would at the same time have $y'''(x)approx -0.05$, so that the second derivative would have to fall, making the first derivative strictly concave and not constant at $y'(x)approx 1$.

In a more dynamical analysis, one would still expect $y'''(x)approx 0$ and $y(x)/xapprox 1$, so that $0approx xy''(x)+0.1$ which solves, in some even more handwaving fashion, to $y'approx -0.1ln x+c$, $y=-0.1x(ln x-1)+cx+d$. Imposing the infinite boundary condition at some largish $T$, this gives for $x>T$
begin{align}
y''(x)&=-frac{0.1}{x}
y'(x)&=1-0.1ln(x/T), \
y&=1.1x-0.1xln(x/T)+y(T)-1.1T
end{align}

The numerical computations below seem to confirm these estimates, especially the one for $y'$ seems to be an upper boundary. While the second and third derivative tend to zero, the first is continuously falling away from $1$.

Also visible, or rather not visible, is any sign of convergence of the solutions for rising $T$. Neither in the graphs of $y$ in the top plot nor in the initial values for the second derivative in the bottom plot.

def flow(a, t): return spi.odeint(lambda u,t:[u[1],u[2],-0.5*u[0]*u[2]-0.05],[0,0,a],t,atol=1e-4,rtol=1e-8,mxstep=2000);



fig, ax = plt.subplots(3,1)

t_plot = np.arange(0,40.01,0.05);

for T in [5,8,12,15,20,25]:

    a = spo.fsolve(lambda x: flow(x,[0,T])[-1,1]-1, 0.1)

    y_plot = flow(a,t_plot);

    for k in range(3): ax[k].plot(t_plot, y_plot[:,k], label="T=%.1f"%T);

for k in range(3): ax[k].legend(); ax[k].grid(); 

plt.show();

Runge-Kutta can be used to solve systems of first-order ODE. Any high-order ODE can be rewritten as a system of first-order ODE's as is written in the article. The equation system written in article, together with the initial conditions for $y$, $f$ and $g$, can directly be plugged into any forwards-integrator, for example from Scipy, by providing it the right-hand-side of the first order ODE system. You will get an array of time-changing values of $y, f, g$ as an answer