First order ODE

ODEs

An ordinary differential equation is an equation of the form

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: Intuition

The Euler method is a fundamental numerical technique for approximating solutions to ordinary differential equations (ODEs). At its core, the method relies on the geometrical concept of a tangent line.

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)$.

Euler method

Consider explicit case where we can write

$$ \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) $$

Euler method

Consider explicit case where we can write

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

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.

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: Intuition

Heun's method is an improvement over Euler's method, aiming to provide a more accurate approximation to the solution of an ordinary differential equation (ODE). It essentially performs a 'correction' step after the initial Euler estimate.

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.

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

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 this plot.

Assignment 3

Assignment 3 will look at these methods for solving ODEs.

Like with Assignment 2 you want to find a way to implement the equations shown today mathematically.

Remember numpy arrays!