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 at 9:18
• Seems like a homework... It is really bad written though... – Alone Programmer Jan 28 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 at 1:54

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[0] = k - (1 - 5/Y[1])*(3+2/Y[1])
dY[1] = Y[0]
dY[2] = 1/Y[1]**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[1],'-', label='r(t)')
plt.plot(t, Y[2],'-', label='phi(t)')
plt.legend(loc='best')
plt.xlabel('T')