Interactive Lecture
Euler Method, Runge-Kutta Methods, Error Analysis, and Stability
by Dr Arnaldo Delli Carri
Many engineering problems involve ordinary differential equations (ODEs) that cannot be solved analytically. Numerical methods provide approximate solutions by discretising the continuous problem into a sequence of algebraic calculations.
Numerical ODE solvers approximate the solution by taking small steps in time, using the differential equation to predict how the solution evolves. The accuracy and stability of these methods depend on the step size and the mathematical properties of the algorithm.
We seek to solve differential equations of the form:
\[\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0\]
Where:
The Euler method is the simplest numerical method, based on a first-order Taylor series approximation.
Algorithm:
\[y_{n+1} = y_n + h \cdot f(t_n, y_n)\]
Where:
Geometric Interpretation: The method follows the tangent line at \((t_n, y_n)\) for a distance \(h\).
Local Truncation Error: \(O(h^2)\)
Global Error: \(O(h)\) - first-order accurate
Euler's method is simple but has poor accuracy and stability. It requires very small step sizes for acceptable results and can be unstable for stiff equations.
An improvement over Euler's method that uses a predictor-corrector approach:
Predictor:
\[\tilde{y}_{n+1} = y_n + h \cdot f(t_n, y_n)\]
Corrector:
\[y_{n+1} = y_n + \frac{h}{2}\left[f(t_n, y_n) + f(t_{n+1}, \tilde{y}_{n+1})\right]\]
Global Error: \(O(h^2)\) - second-order accurate
Runge-Kutta (RK) methods achieve higher accuracy by evaluating the derivative at multiple intermediate points.
Second-Order Runge-Kutta (RK2):
\[k_1 = h \cdot f(t_n, y_n)\]
\[k_2 = h \cdot f\left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right)\]
\[y_{n+1} = y_n + k_2\]
Global Error: \(O(h^2)\)
The most widely used method, offering an excellent balance of accuracy and computational cost:
RK4 Algorithm:
\[k_1 = h \cdot f(t_n, y_n)\]
\[k_2 = h \cdot f\left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right)\]
\[k_3 = h \cdot f\left(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}\right)\]
\[k_4 = h \cdot f(t_n + h, y_n + k_3)\]
\[y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\]
Global Error: \(O(h^4)\) - fourth-order accurate
RK4 evaluates the derivative at four points: the start, two midpoints, and the end of the interval. It then takes a weighted average, giving more weight to the midpoint evaluations.
Knowing where your numerical errors come from — and how large they are — is what makes a solution trustworthy:
Local Truncation Error (LTE): Error introduced in a single step
Global Error: Accumulated error over many steps
Order of Accuracy:
Absolute Error: \(|y_{true} - y_{numerical}|\)
Relative Error: \(\frac{|y_{true} - y_{numerical}|}{|y_{true}|}\)
Stability determines whether errors grow or decay as the solution progresses:
Test Equation: \(\frac{dy}{dt} = \lambda y\) where \(\lambda < 0\)
Stability Criterion: For Euler's method:
\[|1 + h\lambda| \leq 1 \implies h \leq \frac{2}{|\lambda|}\]
Stiff Equations: Equations with widely different time scales (large negative \(\lambda\)) require very small step sizes for explicit methods or implicit methods for stability.
Modern ODE solvers adjust step size automatically to maintain accuracy:
Embedded RK Methods: Use two approximations of different orders
Error estimate: \(\epsilon = |y_{higher} - y_{lower}|\)
If \(\epsilon > \text{tolerance}\): reduce \(h\)
If \(\epsilon < \text{tolerance}\): increase \(h\)
All methods extend naturally to systems of equations:
\[\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y}), \quad \mathbf{y}(t_0) = \mathbf{y}_0\]
where \(\mathbf{y}\) and \(\mathbf{f}\) are vectors. Apply the method component-wise.
Example: Second-order ODE \(\frac{d^2y}{dt^2} = g(t, y, \frac{dy}{dt})\)
Convert to system: \(y_1 = y, \quad y_2 = \frac{dy}{dt}\)
\[\frac{dy_1}{dt} = y_2, \quad \frac{dy_2}{dt} = g(t, y_1, y_2)\]
Compare different numerical methods solving the same differential equation. Observe how accuracy varies with step size and method choice.
The absolute error of each method at the final time is overlaid on the plot — read accuracy against cost.
| Method | Algorithm | Order | Evaluations/Step | Stability |
|---|---|---|---|---|
| Euler | \(y_{n+1} = y_n + hf(t_n, y_n)\) | 1st | 1 | Conditionally stable |
| Modified Euler | \(y_{n+1} = y_n + \frac{h}{2}[f(t_n, y_n) + f(t_{n+1}, \tilde{y}_{n+1})]\) | 2nd | 2 | Better than Euler |
| RK2 | Midpoint method | 2nd | 2 | Good |
| RK4 | \(y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)\) | 4th | 4 | Excellent |
Use Euler's method with step size \(h = 0.2\) to solve \(\frac{dy}{dt} = -2y\) with \(y(0) = 1\) from \(t = 0\) to \(t = 0.4\). Compare with the exact solution \(y(t) = e^{-2t}\).
Show Solution
Given:
Euler's Method: \(y_{n+1} = y_n + h \cdot f(t_n, y_n) = y_n + 0.2(-2y_n) = y_n(1 - 0.4)\)
Step 1: \(t_1 = 0.2\)
\[y_1 = y_0(1 - 0.4) = 1(0.6) = 0.600\]
Exact: \(y(0.2) = e^{-0.4} = 0.6703\)
Error: \(|0.6703 - 0.600| = 0.0703\)
Step 2: \(t_2 = 0.4\)
\[y_2 = y_1(1 - 0.4) = 0.6(0.6) = 0.360\]
Exact: \(y(0.4) = e^{-0.8} = 0.4493\)
Error: \(|0.4493 - 0.360| = 0.0893\)
Summary Table:
| \(t\) | Euler \(y_n\) | Exact | Error |
|---|---|---|---|
| 0.0 | 1.0000 | 1.0000 | 0.0000 |
| 0.2 | 0.6000 | 0.6703 | 0.0703 |
| 0.4 | 0.3600 | 0.4493 | 0.0893 |
Error grows with time — smaller step sizes would improve accuracy.
Answer: \(y(0.2) = 0.600\) (exact \(0.6703\), error \(0.0703\)); \(y(0.4) = 0.360\) (exact \(0.4493\), error \(0.0893\)).
Apply one step of RK4 with \(h = 0.5\) to solve \(\frac{dy}{dt} = t + y\) with \(y(0) = 1\).
Show Solution
Given:
Calculate \(k_1\):
\[k_1 = h \cdot f(t_0, y_0) = 0.5(0 + 1) = 0.5\]
Calculate \(k_2\):
\[k_2 = h \cdot f\left(t_0 + \frac{h}{2}, y_0 + \frac{k_1}{2}\right) = 0.5 \cdot f(0.25, 1.25)\]
\[k_2 = 0.5(0.25 + 1.25) = 0.5(1.5) = 0.75\]
Calculate \(k_3\):
\[k_3 = h \cdot f\left(t_0 + \frac{h}{2}, y_0 + \frac{k_2}{2}\right) = 0.5 \cdot f(0.25, 1.375)\]
\[k_3 = 0.5(0.25 + 1.375) = 0.5(1.625) = 0.8125\]
Calculate \(k_4\):
\[k_4 = h \cdot f(t_0 + h, y_0 + k_3) = 0.5 \cdot f(0.5, 1.8125)\]
\[k_4 = 0.5(0.5 + 1.8125) = 0.5(2.3125) = 1.15625\]
Final Update:
\[y_1 = y_0 + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\]
\[y_1 = 1 + \frac{1}{6}(0.5 + 2(0.75) + 2(0.8125) + 1.15625)\]
\[y_1 = 1 + \frac{1}{6}(0.5 + 1.5 + 1.625 + 1.15625)\]
\[y_1 = 1 + \frac{4.78125}{6} = 1 + 0.7969 = 1.7969\]
The exact solution is \(y = 2e^t - t - 1\), giving \(y(0.5) = 1.7974\), so the absolute error is \(\approx 0.0006\).
Answer: \(y(0.5) \approx 1.797\) (absolute error \(\approx 0.0006\)).
For the equation \(\frac{dy}{dt} = -10y\) with \(y(0) = 1\):
(a) What is the maximum step size for Euler's method to be stable?
(b) Compute \(y(0.2)\) using step sizes \(h = 0.05\) and \(h = 0.025\), and verify that the error ratio is consistent with Euler being a first-order method.
Show Solution
Given:
(a) Stability Condition:
For Euler's method to be stable: \(|1 + h\lambda| \leq 1\)
\[|1 + h(-10)| \leq 1\]
\[|1 - 10h| \leq 1\]
This gives: \(-1 \leq 1 - 10h \leq 1\)
From the left inequality: \(-1 \leq 1 - 10h \implies 10h \leq 2 \implies h \leq 0.2\)
From the right inequality: \(1 - 10h \leq 1\) (always satisfied for \(h > 0\))
So the maximum stable step size is \(h_{max} = 0.2\).
(b) Step-Size Convergence Study:
For Euler's method, the amplification factor per step is \(R = 1 + h\lambda = 1 - 10h\).
After \(n\) steps to reach \(t = 0.2\): \(y_n = R^n\).
With \(h = 0.05\) (4 steps): \(R = 1 - 0.5 = 0.5\)
\[y_4 = (0.5)^4 = 0.0625\]
Exact: \(y(0.2) = e^{-2} = 0.1353\)
Error: \(|0.1353 - 0.0625| = 0.0728\)
With \(h = 0.025\) (8 steps): \(R = 1 - 0.25 = 0.75\)
\[y_8 = (0.75)^8 = 0.1001\]
Error: \(|0.1353 - 0.1001| = 0.0352\)
Error ratio: \(\frac{0.0728}{0.0352} \approx 2.07 \approx 2\)
Halving the step size roughly halved the error, confirming Euler is first-order accurate (\(O(h)\)).
For an \(n\)th-order method, halving \(h\) reduces the global error by a factor of \(2^n\). This step-size halving test is a practical way to verify the order of any numerical method.
Answer: (a) \(h_{max} = 0.2\); (b) error ratio \(\approx 2.07 \approx 2\), confirming Euler is first-order accurate (\(O(h)\)).
Use RK4 with \(h = 0.1\) to take one step for the system:
\[\frac{dx}{dt} = y, \quad \frac{dy}{dt} = -x, \quad x(0) = 1, \quad y(0) = 0\]
Show Solution
Given:
Compute all \(k\) values for both components simultaneously:
\(k_1\):
\[k_1^x = h \cdot f_1(x_0, y_0) = 0.1(0) = 0\]
\[k_1^y = h \cdot f_2(x_0, y_0) = 0.1(-1) = -0.1\]
\(k_2\): Evaluate at \((x_0 + k_1^x/2,\; y_0 + k_1^y/2) = (1,\; -0.05)\)
\[k_2^x = h \cdot f_1(1, -0.05) = 0.1(-0.05) = -0.005\]
\[k_2^y = h \cdot f_2(1, -0.05) = 0.1(-1) = -0.1\]
\(k_3\): Evaluate at \((x_0 + k_2^x/2,\; y_0 + k_2^y/2) = (0.9975,\; -0.05)\)
\[k_3^x = h \cdot f_1(0.9975, -0.05) = 0.1(-0.05) = -0.005\]
\[k_3^y = h \cdot f_2(0.9975, -0.05) = 0.1(-0.9975) = -0.09975\]
\(k_4\): Evaluate at \((x_0 + k_3^x,\; y_0 + k_3^y) = (0.995,\; -0.09975)\)
\[k_4^x = h \cdot f_1(0.995, -0.09975) = 0.1(-0.09975) = -0.009975\]
\[k_4^y = h \cdot f_2(0.995, -0.09975) = 0.1(-0.995) = -0.0995\]
Update:
\[x_1 = 1 + \frac{1}{6}(0 + 2(-0.005) + 2(-0.005) + (-0.009975)) = 1 - 0.004996 = 0.99500\]
\[y_1 = 0 + \frac{1}{6}(-0.1 + 2(-0.1) + 2(-0.09975) + (-0.0995)) = -0.09983\]
Exact Solution:
\[x(0.1) = \cos(0.1) = 0.99500\]
\[y(0.1) = -\sin(0.1) = -0.09983\]
Comparing against the exact values, the absolute errors after one step are \(|x_1 - \cos(0.1)| \approx 1.4 \times 10^{-9}\) and \(|y_1 + \sin(0.1)| \approx 8.3 \times 10^{-8}\) — RK4 is accurate to roughly eight significant figures here.
Answer: \(x_1 = 0.99500\), \(y_1 = -0.09983\); absolute errors \(\approx 1.4 \times 10^{-9}\) and \(\approx 8.3 \times 10^{-8}\).