Solving ODEs Numerically

Many physical systems are described by differential equations that cannot be solved analytically. This is common in systems with many interacting components, such as gravitational N-body problems or fluid dynamics in weather and climate models. In these cases, we turn to numerical methods to approximate the solution.

Suppose we have a collection of time-dependent variables $ y_i(t) $, which may represent quantities such as position, velocity, temperature, or electric field. If we can express how each variable evolves in time, we can write a system of first-order differential equations:

$$ \frac{dy_0}{dt} = F_0(y_0, y_1, \dots, y_{n-1}, t), $$ $$ \frac{dy_1}{dt} = F_1(y_0, y_1, \dots, y_{n-1}, t), $$ $$ \dots $$ $$ \frac{dy_{n-1}}{dt} = F_{n-1}(y_0, y_1, \dots, y_{n-1}, t) $$

The functions $F_i$ describe how the system evolves and may depend on both the variables themselves and time. This formulation is general and applies to a wide range of physical systems.

Numerical methods such as Euler's method or Runge-Kutta schemes approximate the solution by advancing the system in small time steps. Starting from initial conditions at time $t_k$, a simple update step takes the form

$$ y_i(t_{k+1}) = y_i(t_k) + \Delta t \, F_i(y_0(t_k), \dots, y_{n-1}(t_k)). $$

The derivative determines the local rate of change, so numerical methods repeatedly use this information to estimate the future state of the system.

In practice, modern libraries implement more sophisticated versions of these methods. In Python, the SciPy function solve_ivp() provides a flexible interface for solving systems of differential equations.

To use solve_ivp(), the system is written as a function that returns all derivatives for a given state and time. The state of the system is passed as an array containing the values of the variables:


          def dydt(t, y):
              # Compute derivatives at time t
              return [dy0_dt, dy1_dt, dy2_dt, ...]

          # Initial values for variables y0, y1, y2, ... 
          y0 = [y0_0, y1_0, ..., yn_0]

          # The time steps at which you want to evaluate the solution
          t_eval = np.linspace(t_0, t_final, N)

          sol = solve_ivp(
              dydt, (t_0, t_final), y0,
              t_eval = t_eval,
              method = "LSODA"
          )

The solution object contains the time points and corresponding values of the variables:


          t_sol = sol["t"]
          y_sol = sol["y"]

Here, y_sol is a two-dimensional array. Each row corresponds to one variable, while each column corresponds to a time value.

The argument method = "LSODA" selects an adaptive algorithm that automatically switches between stiff and non-stiff methods, making it a reliable general-purpose choice.

Two example systems are shown below. Both are non-dimensionalized equations that we will solve using solve_ivp()

Example 1: Non-dimensionalized free fall with drag (see derivation)

$$ \frac{du}{d\tau} = 1 - u^2 $$

Example 2: Non-dimensionalized damped harmonic oscillator (see derivation)

$$ \frac{d^2 l}{d\tau^2} + \gamma \frac{dl}{d\tau} + l = 0 $$ Open Notebook Environment