5
$\begingroup$

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

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Federico Poloni Jan 28 at 9:18
  • $\begingroup$ Seems like a homework... It is really bad written though... $\endgroup$ – Alone Programmer Jan 28 at 19:21
  • $\begingroup$ @Лада, I've edited the formulas in your question (btw, we have LaTeX enabled!). Please, make sure that I did not miss anything. $\endgroup$ – Anton Menshov Jan 29 at 1:54
2
$\begingroup$

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')
| cite | improve this answer | |
$\endgroup$
6
$\begingroup$

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

resulting plot

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ nice answer but I think that the author has mistaken the square of a first-order derivative for a second-order derivative $\endgroup$ – QuantumApple Apr 27 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.