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

  • 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

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

resulting plot


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.