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$$
$$ 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.
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)$.
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) $$
Consider explicit case where we can write
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} $$
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.
The algorithm is:
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 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 $$
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 this plot.
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!