First order ODE

ODEs

An ordinary differential equation is an equation of the form

$$G(t,x(t),x’(t),x’’(t),…)=0$$

  • It is called first order if only $x$ and $x’$ appear in $G$.
  • a first order ODE is called explicit if one can write it as

$$ x’(t) = f(x(t),t)$$

Differential equation applications

Differential equations describe the dynamics of many systems:

  • Physics:
    • ballistics, planets and spacecrafts trajectories (Newton’s law)
    • fluid dynamics (Navier-Stokes equation)
    • Electromagnetism (Maxwell’s equations)
    • General relativity (Einstein’s field equation)
    • Thermodynamics
  • Chemistry
  • Biology
  • Economics

In most cases the equations cannot be solved analytically: we need numerical methods.

ODE vs integration

Integration and ODEs are related.

Integration: $$ \frac{dx}{dt} = f(t)$$

Differential equation: $$ \frac{dx}{dt} = f(x(t),t)$$

An ODE where $f$ does not depend on $x(t)$ is an integration.

Euler method

Key Concept:

If you know the current state $x(t)$ and the rate of change $f(x(t), t)$ at that state, you can approximate the state at the next time step $t + \Delta t$ by 'walking' along the tangent line.

This approximation assumes that $f(x(t), t)$ remains relatively constant over a small time step $\Delta t$, allowing us to linearly extrapolate $x(t)$ to $x(t + \Delta t)$.

$$ \frac{dx}{dt}=f(x(t),t)$$

We start from the usual place: the definition of the derivative.

$$ \frac{dx}{dt} = \lim_{\Delta t\rightarrow 0} \frac{x(t+\Delta t)-x(t)}{\Delta t}$$

Solve for $x(t+\Delta t)$:

$$ x(t+\Delta t) \approx x(t) + \Delta t \frac{dx}{dt} = x(t) + \Delta t f(x(t),t) $$

Algorithmic Steps:

  1. Initialize $x_0$ and $t_0$ as initial conditions.
  2. Choose a time step $\Delta t$.

Iterative Update:

$$ x_{n+1} = x_n + \Delta t \cdot f(x_n, t_n) $$ $$ t_{n+1} = t_n + \Delta t $$

Example: Radioactive decay

  • A population of unstable nuclei undergoes radioactive decay
  • A single nuclei decays at a random time
  • The ‘continuum behaviour’ of this random process can be described by an ODE
    • continuum means we treat the number of nuclei as a real number instead of an integer
    • statistically meaningful but not realistic for a single experiment

Number of nuclei at time $t_0$: $N_0$

Mean lifetime of the decay process: $\tau$

Differential equation:

$$ \frac{dN}{dt} = - \frac{N(t)}{\tau}\equiv f(N,t)$$

Analytical solution:

$$ N(t)= N_0 e^{-t/\tau} $$

In Assignment 3, be careful not to use your analytical solution in the differential equation $f_{rad}$ which is a function of $N$ and $t$: $N_0 \neq N \quad \forall t$!

With the analytical solution we can relate the derivative and the number of nuclei:

The area under the curve is the change in the number of nuclei.

If we do not know the analytical solution things get more complicated, we need to:

  • know $N$ to calculate the derivative
  • know the derivative to calculate the next position

So we get the algorithm

  • estimate the area under the derivative for a small step
  • change position according to this area
  • repeat

Different methods have different ways of estimating the area under the derivative curve. Euler’s method estimates the area under the derivative using the rectangle method

order integration method ODE method
0 Rectangle Euler
1 Trapezium rule Heun
2 Simpson rule Runge Kutta 4

Euler’s method demo

Heun’s method

Heun’s method uses the same idea than the trapezium method to improve the estimate of the area under the derivative curve.

Key Concept:

Instead of using just the initial slope to advance to the next point, Heun's method also calculates the slope at the new point. It then takes the average of these two slopes to make a more accurate step forward, mimicking the area under the curve more closely.

The algorithm is:

  • use the slope at the current point to jump to the next point
  • calculate the slope there
  • use the average of the two slopes to advance

Heun's Method: Mathematics

Heun's method is formally defined as follows:

Let $ \frac{dx}{dt} = f(x(t), t) $. The update equation is:

$$ x(t+\Delta t) \approx x(t) + \frac{\Delta t}{2} [ f(x(t), t) + f(x(t) + \Delta t f(x(t), t), t+\Delta t) ] $$

Heun's Method: Algorithm

Algorithmic Steps:

  1. Initialize $ x_0 $ and $ t_0 $ as initial conditions.
  2. Choose a time step $ \Delta t $.

Iterative Update:

$$ x_{n+1} = x_n + \frac{\Delta t}{2} [ f(x_n, t_n) + f(x_n + \Delta t f(x_n, t_n), t_n+\Delta t) ] $$ $$ t_{n+1} = t_n + \Delta t $$

Heun’s method demo

Runge Kutta

Heun’s method is one of the Runge Kutta family of methods that use multiple intermediate steps to determine the value of the function after a time step.

The most used is one involving four intermediate steps, it is called Runge-Kutta 4 but is it often simply named “Runge-Kutta”

Intuition Behind Runge-Kutta Method

The essence of RK4 lies in its strategic choice of multiple sampling points within each time step. Unlike simpler methods, RK4 uses these points to compute a weighted average of derivative estimates.

Why It Matters: This approach enables RK4 to capture more complex behaviors of the function, achieving higher accuracy without drastically increasing computational cost.

Mathematical Overview of Runge-Kutta

Given an ODE of the form $$ \frac{dx}{dt} = f(x(t), t), $$ the RK4 method calculates the next step as follows:

$$ x(t + \Delta t) = x(t) + \frac{\Delta t}{6} \left(k_1 + 2k_2 + 2k_3 + k_4 \right) $$

Algorithmic Steps for Runge-Kutta

  1. Initialize initial conditions $$ x_0, t_0 $$
  2. Choose time step $$ \Delta t $$
  3. Iterate using: $$ k_1 = f(x_n, t_n) ,$$ $$ k_2 = f(x_n + \frac{\Delta t}{2} k_1, t_n + \frac{\Delta t}{2}) ,$$ $$ k_3 = f(x_n + \frac{\Delta t}{2} k_2, t_n + \frac{\Delta t}{2}) ,$$ $$ k_4 = f(x_n + \Delta t k_3, t_n + \Delta t) ,$$ $$ x_{n+1} = x_n + \frac{\Delta t}{6} (k_1 + 2k_2 + 2k_3 + k_4) ,$$ $$ t_{n+1} = t_n + \Delta t .$$

The steps are :

  • calculate the slope at the start, call it $k_1$
  • use $k_1$ to advance from the start to the middle of the time step interval
  • calculate the slope there, call it $k_2$
  • use $k_2$ to advance from the start to the middle of the time step interval
  • calculate the slope there, call it $k_3$
  • use $k_3$ to advance from the start to the end of the time step interval
  • calculate the slope there, call it $k_4$
  • combine the four slopes as

$$ k = \frac{1}{6}\left(k_1+2k_2+2 k_3 + k_4\right)$$

  • use $k$ as the slope to advance to the next step

Runge Kutta 4 demo

Error evaluation

ODE method # function evaluations error per step total error
Euler 1 $\mathcal{O}(\Delta t)^2$ $\mathcal{O}(\Delta t)^1$
Heun 2 $\mathcal{O}(\Delta t)^3$ $\mathcal{O}(\Delta t)^2$
Runge Kutta 4 4 $\mathcal{O}(\Delta t)^5$ $\mathcal{O}(\Delta t)^4$

The error scaling is shown in the plot below; here we look at the difference at 10,000 years. For Assignment 3 your error plot does not have to be identical to this. Don't worry about slight differences.

np.logspace has arguments equivalent to (starting power, ending power, number of values) this means if we want to go from $10^{0}$ to $10^{4}$ in 5 steps we can write "np.logspace(0, 4, 5, dtype=int)" which will create a numpy array $[1, 10, 100, 1000, 10000]$.