Solving coupled differential equations in Python, 2nd order

Question

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} $$

Answer 1

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[0] 
    f_r = y[1] # 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')

Answer 2

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')