# Solving coupled differential equations in Python, 2nd order

I have a system of coupled differential equations, one of which is second-order. I am looking for a way to solve them in Python. I would be extremely grateful for any advice on how can I do that!

$$k$$ is just a constant

$$\left(\frac{dr}{dT}\right)^2=k-\left(1-\frac{5}{r}\right)\left(3+\frac{2}{r^2}\right)\\ \frac{d\varphi}{dT}=\frac{1}{r^2}$$

• Is that $(\frac{dr}{dt})^2$ the square of $\frac{dr}{dt}$, or a second derivative $\frac{d^2}{dt^2}r$? That notation is typically used for the first derivative squared, but you mention "second order" twice, so it is not clear. – Federico Poloni Jan 28 '20 at 9:18
• Seems like a homework... It is really bad written though... – Alone Programmer Jan 28 '20 at 19:21
• @Лада, I've edited the formulas in your question (btw, we have LaTeX enabled!). Please, make sure that I did not miss anything. – Anton Menshov Jan 29 '20 at 1:54

This is a system of first order differential equations, not second order. It models the geodesics in Schwarzchield geometry. In other words, this system represents the general relativistic motion of a test particle in static spherically symmetric gravitational field. In general, there is a third equation for how coordinate time is related to proper time.

The full system should be \begin{align} &\left(\frac{dr}{dT}\right)^2 =\frac{E^2}{m^2c^2} \, - \, \left(1 - \frac{r_{s}}{r}\right)\left(c^2 + \frac{h^2}{r^2}\right)\\ &\\ &\frac{d\varphi}{dT} = \frac{h}{r^2}\\ &\\ &\frac{dt}{dT} = \frac{E}{mc^2}\left(\frac{r}{r - r_s}\right) \end{align} where $$m$$ is the mass of the test particle, $$E$$ is the energy of the particle, $$r_s$$ is the Schwarzschield radius, and $$c$$ is the speed of light in vacuum. I do not know why in your case you have $$h = \sqrt{2}$$ from the first equation but $$h = 1$$ from the second, so you should check that. Either way, your system looks like this \begin{align} &\left(\frac{dr}{dT}\right)^2 = k \, - \, \left(1 - \frac{5}{r}\right)\left(3 + \frac{2}{r^2}\right)\\ &\\ &\frac{d\varphi}{dT} = \frac{1}{r^2}\\ &\\ &\frac{dt}{dT} = k_1\left(\frac{r}{r - 5}\right) \end{align}

For a quick and dirty calculation, I would differentiate the first equation once with respect to $$T$$: \begin{align} &\frac{d}{dT}\left(\frac{dr}{dT}\right)^2 = \frac{d}{dT}\left( \, k \, - \, \left(1 - \frac{5}{r}\right)\left(3 + \frac{2}{r^2}\right) \, \right)\\ &2 \, \frac{dr}{dT} \, \frac{d^2r}{dT^2} = -\, \frac{d}{dT}\left( \, \left(1 - \frac{5}{r}\right)\left(3 + \frac{2}{r^2}\right) \, \right)\\ &2 \, \frac{d^2r}{dT^2} \, \frac{dr}{dT} = -\, \frac{d}{dT}\left( \, 3 - \frac{15}{r} + \frac{2}{r^2} - \frac{10}{r^3}\, \right)\\ &2 \, \frac{d^2r}{dT^2} \, \frac{dr}{dT} = -\, \left( \, + \frac{15}{r^2} - 2\, \frac{2}{r^3} + 3\, \frac{10}{r^4}\, \right) \frac{dr}{dT} \end{align} Cancel out $$\frac{dr}{dT}$$ on both sides of the equation and open the parentheses: \begin{align} &2 \, \frac{d^2r}{dT^2} = - \frac{15}{r^2} + \frac{4}{r^3} - \frac{30}{r^4} \end{align} Now this is a second order differential equation. So if you introduce the variable $$u = \frac{dr}{dT}$$ you get the system \begin{align} &\frac{dr}{dT} = u\\ &\\ &\frac{du}{dT} = - \frac{7.5}{r^2} + \frac{2}{r^3} - \frac{15}{r^4} \\ &\\ &\frac{d\varphi}{dT} = \frac{1}{r^2}\\ &\\ &\frac{dt}{dT} = k_1\left(\frac{r}{r - 5}\right) \end{align} Given some initial conditions $$r_0, \, \varphi_0, \, t_0$$ for the original system, you have to calculate an initial condition for the new variable $$u$$. To that end, you have to calculate the equation: $$u_0 = \pm \sqrt{k \, - \, \left(1 - \frac{5}{r_0}\right)\left(3 + \frac{2}{r_0^2}\right)}$$

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

# y = np.array([r, u, phi, time])
def f(t, y):
r = y
f_r = y # this is the dr / dT auxiliary equation
f_u = - 7.5 / (r**2) + 2 / (r**3) - 15 / (r**4)
f_phi = 1 / (r**2)
f_time = k1 * r / (r - 5) # this is the equation of the time coordinate
return np.array([f_r, f_u, f_phi, f_time])

# from the initial value for r = r0 and given energy k,
# calculate the initial rate of change dr / dT = u0
def ivp(r0, k, sign):
u0 = math.sqrt( k - ( 1 - 5 / (r0**2) ) * ( 3 + 2 / (r0**2) ) )
return sign * u0

k = 3.0
k1 = 2.0
r0 = 20.0
sign = 1 # or -1

u0 = ivp(r0, k, sign)
# y = np.array([r, u, phi, time])
y0 = [r0, u0, math.pi/6, 0]
t_span = np.linspace(0, 1000, num=1001)

sol = solve_ivp(f, [0, 1000], y0, method='Radau', t_eval=t_span)

plt.plot(sol.t, sol.y[0,:],'-', label='r(t)')
plt.plot(sol.t, sol.y[2,:],'-', label='phi(t)')
plt.legend(loc='best')
plt.xlabel('T')


The first step is to transform the second order equation to a set of two coupled first order equations. Define an auxiliary function $$u(T) = \frac{dr(T)}{dT}$$. This results in the system

\begin{align} \frac{du}{dT} &= k-(1-\frac{5}{r})(3+\frac{2}{r^2}) \\ \frac{dr}{dT} &= u\\ \frac{d\phi}{dT} & = \frac{1}{r^2} \end{align}

Now you have a set of three coupled first order equations in the form fit for solving with solve_ivp. See SciPy documentation for solve_ivp.

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

def rhs(t,Y):
dY = np.zeros_like(Y)
k = 1.0
dY = k - (1 - 5/Y)*(3+2/Y)
dY = Y
dY = 1/Y**2
return dY

Y0 = np.array([0,1,0])
sol = solve_ivp(rhs, [0,10], Y0, method='Radau', dense_output=True)

t = np.linspace(0, 10, 1001)
Y = sol.sol(t)

plt.plot(t, Y,'-', label='r(t)')
plt.plot(t, Y,'-', label='phi(t)')
plt.legend(loc='best')
plt.xlabel('T') • nice answer but I think that the author has mistaken the square of a first-order derivative for a second-order derivative – QuantumApple Apr 27 '20 at 20:19