# Solving second order ODE numerically

I have the following second order differential equation I want to solve numerically in Python (or Matlab):

$$$$\frac{d^2y}{dx^2}=a \left[ \left(\frac{y}{b}\right)^{-3} - \left(\frac{y}{b}\right)^{-6} \right]$$$$ with initials conditions $$y(0)=b$$ and $$\frac{dy}{dx}(0)=c$$, where where $$a$$,$$b$$,$$c$$ are some constants.

Now I reduced it to 2 first order ODEs when setting $$p_1=\frac{dy}{dx}$$ and $$p_2=y(x)$$:

$$$$\frac{dp_1}{dx}=a \left[ \left(\frac{y}{b}\right)^{-3} - \left(\frac{y}{b}\right)^{-6} \right], \hspace{3em} \frac{dp_2}{dx}=p_1$$$$

But how to continue from here?

• Are you sure this is a scalar second order DE and not the field strength of a central field, similar to the Lennard-Jones potential? Aug 9 at 11:54
• @LutzLehmann. It really look like the LJ potential. Good catch. Cheers :-) Aug 9 at 12:16

Here is how I solved it in MATLAB:

h0 = 1;  % initial value for h
dhdx0 = 1e-5;  % initial value for dh/dx

xspan = [0 1e4];
y0 = [h0 dhdx0];

c = (5/3) * (1e-3)^2;

[x,y] = ode45(@(x,y) odefcn(x,y,c), xspan, y0);

plot(x,y(:,1),'-o',x,y(:,2),'-.')

function dydt = odefcn(x,y,c)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = c * (y(1)^(-3) - y(1)^(-6));
end


I tried to do the same in Python, but get strange results (think that MATLABs smart mesh is the difference):

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

h0 = 1  # initial value for h
dhdx0 = 1e-5  # initial value for dh/dx

m = GEKKO(remote=False)
m.time = np.linspace(0, 1e4, 250)
y = m.Var(h0)
dy = m.Var(dhdx0)

c = (5/3) * (1e-3)**2
m.Equation(y.dt() == dy)
m.Equation(dy.dt() == c * (y**(-3) - y**(-6)))

m.options.IMODE = 6
m.options.NODES = 3
m.solve(disp=False)

plt.plot(m.time, y.value, '.-')
plt.show()


EDIT: as suggested by @LutzLehmann, the equivalent to ODE45 in Python is scipy.integrate.solve_ivp:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

h0 = 1  # initial value for h
dhdx0 = 1e-5  # initial value for dh/dx

def f(x, y, c):
return y[1], c * (y[0]**(-3) - y[0]**(-6))

xspan = np.linspace(0, 1e4, 100)
c = (5/3) * (1e-3)**2
sol = solve_ivp(lambda x, y: f(x, y, c), [xspan[0], xspan[-1]], [h0, dhdx0], t_eval=xspan, rtol=1e-5)

plt.plot(sol.t, sol.y[0], '.-', label='Output')
plt.show()

• The direct equivalent in python is scipy.integrate.solve_ivp with the standard method RK45, there you also get the adaptive mesh. For differences take notice of the Refine option in matlab's solvers that does not exist in scipy. Aug 12 at 6:36
• @LutzLehmann, thanks! That solved the problem. Solution updated. Aug 13 at 9:59