Numerical methods for solving coupled ODEs?

I am currently learning differential geometry and general relativity on my own. I have been running through exercises to practice calculating geodesics. However, I notice that most of the ODEs from the geodesic equation are impossible to solve analytically. Which numerical methods work to solve these types of equations (example below)?

For example, the equations for hyperbolic space are:

$\frac{d^{2}u}{d\lambda^{2}}+\sinh(2u)(\frac{dv}{d\lambda})^{2}=0$

$\frac{d^{2}v}{d\lambda^{2}}+2\coth(u)\frac{du}{d\lambda}\frac{dv}{d\lambda}=0$

Again, in summary, what numerical methods work to solve equations like these?

Solution 1:

Google "Rewrite higher order ODE as system" for some guides (e.g. This one. There's also some examples in the matlab docs for their ode45 ODE solver.).

Basically we just use the fact that $\mathrm{D}^{n+1} f = \mathrm{D} ( \mathrm{D}^n f)$ and define all the intermediate steps as new functions that we then connect via the actual equation. Let's say you have an explicit ODE like $$\mathrm{D}^n u (t) = f(t, u, \mathrm{D} u, ..., \mathrm{D}^{n-1} u),$$ then you can turn it into the equivalent system $$\begin{align*} \mathrm{D} x_i &= x_{i+1} \quad i=1,...,n-1\\ \mathrm{D} x_n(t) &= f(t, x_1, ..., x_n) \end{align*}$$ where a solution $(x_1, ..., x_n)$ satisfies $x_1 = u, x_2 = \mathrm{D} u, ..., x_n = \mathrm{D}^n u$.

So for example: consider the geodesics on the hyperbolic paraboloid parametrized by $$\phi: U \to \Bbb R^3, (u,v) \mapsto \phi(u,v):=(u, v, \frac{u^2-v^2}{2})^T.$$ The geodesic equation for this case simplifies to

$$\begin{align*} \mathrm{D}^2 y^u + \frac{u}{u^2+v^2+1}((\mathrm{D} y^u)^2 - (\mathrm{D} y^v)^2) &= 0 \\ \mathrm{D}^2 y^v + \frac{v}{u^2+v^2+1}((\mathrm{D} y^v)^2 - (\mathrm{D} y^u)^2) &= 0 \end{align*}$$ where $y=(y^u, y^v)^T$ is a curve on the parameter space $U \subseteq \Bbb R^2$ and $u,v$ are the points through which this curve passes, so precisely $y^u, y^v$. We rewrite this as the following system: $$\begin{align*} \mathrm{D} x_{1,u} &= x_{2,u} \\ \mathrm{D} x_{1,v} &= x_{2,v} \\ \mathrm{D} x_{2,u} &= -\frac{x_{1,u}}{x_{1,u}^2+x_{1,v}^2+1}(x_{2,u}^2 - x_{2,v}^2) \\ \mathrm{D} x_{2,v} &= -\frac{x_{1,v}}{x_{1,u}^2+x_{1,v}^2+1}(x_{2,v}^2 - x_{2,u}^2). \end{align*}$$

Depending on the language you want to use, solving this could for example look like this (this uses my own ODE python library https://github.com/SV-97/odespy/ - but it'll be similar in other languages / with other libraries):

import odes as ode
import numpy as np

def phi3(u, v):
    return (u**2 - v**2)/2

alpha = 0
x_0 = np.array([1, 1, np.cos(alpha), np.sin(alpha)]) # start position is (1, 1) and direction is given by alpha

@ode.solve_ivp(solver=ode.explicit_rk4, t_end=40.0, step_size=1e-2)
@ode.ivp(t_0, x_0)
def geodesic_equation(t: np.array, y: np.array):
    y_u, y_v, dy_u, dy_v = y
    r = y_u**2 + y_v**2 + 1
    h = dy_u**2 - dy_v**2
    d2y_u = -y_u / r * h
    d2y_v = y_v / r * h
    return np.array([dy_u, dy_v, d2y_u, d2y_v])
    
sol = geodesic_equation()
u = sol.ys[:, 0]
v = sol.ys[:, 1]

# now compute the actual geodesic on the surface
x = u
y = v
z = phi3(u, v)