An ordinary differential equation is an equation of the form
$$G(t,x(t),x’(t),x’’(t),…)=0$$
$$ x’(t) = f(x(t),t)$$
Differential equations describe the dynamics of many systems:
In most cases the equations cannot be solved analytically: we need numerical methods.
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.
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:
Iterative Update:
$$ x_{n+1} = x_n + \Delta t \cdot f(x_n, t_n) $$ $$ t_{n+1} = t_n + \Delta t $$
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:
So we get the algorithm
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 |
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:
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) ] $$
Algorithmic Steps:
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 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”
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.
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) $$
The steps are :
$$ k = \frac{1}{6}\left(k_1+2k_2+2 k_3 + k_4\right)$$
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]$.