How to numerically solve this ODE: $y''' \left(x\right)=-0.5\cdot y \left(x\right) \cdot y''\left( x \right) -0.05 $

Question

I am trying to numerically solve the following ordinary differential equation that I encountered in one article: $$ y''' \left( x \right) =-0.5\cdot 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.5\cdot 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.

Answer 1

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

Answer 2

An equilibrium of the demanded type seems improbable due to the constant $-0.05$. If indeed for $x\approx\infty$ 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)/x\approx 1$, so that $0\approx xy''(x)+0.1$ which solves, in some even more handwaving fashion, to $y'\approx -0.1\ln 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.1\ln(x/T), \\ y&=1.1x-0.1x\ln(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();