Euler's Method and Runge–Kutta (RK4) with Worked Examples
When a Differential Equation Has No Neat Solution
Many differential equations that model real systems, such as the swing of a large pendulum, the spread of a disease or the flight of a rocket, cannot be solved with a formula. Engineers and scientists solve them numerically instead: starting from a known point, they take small steps and estimate how the solution changes. The simplest step-by-step method is Euler's method; the most widely used is the fourth-order Runge–Kutta method (RK4). This tutorial explains both, works through them by hand, compares their accuracy against an exact solution, and shows a short Python program that does the arithmetic. It suits MTH 2xx numerical analysis courses and engineering mathematics.
The Problem We Want to Solve
An initial value problem gives a first-order differential equation and a starting point:
dy/dx = f(x, y), with y = y₀ when x = x₀.
We want approximate values of y at x₀ + h, x₀ + 2h, x₀ + 3h and so on, where h is the step size.
Throughout this tutorial we use the test problem dy/dx = x + y, y(0) = 1. It can be solved exactly (it is a linear equation, solvable with an integrating factor): y = 2eˣ − x − 1. Having the exact answer lets us measure the error of each method.
Euler's Method
At the current point (xₙ, yₙ), the gradient of the solution curve is f(xₙ, yₙ). Euler's method follows the tangent line for one step:
yₙ₊₁ = yₙ + h f(xₙ, yₙ), and xₙ₊₁ = xₙ + h.
Example 1. Use Euler's method with h = 0.1 to estimate y(0.2) for dy/dx = x + y, y(0) = 1.
- Step 1: f(0, 1) = 0 + 1 = 1. y₁ = 1 + 0.1 × 1 = 1.1 at x = 0.1.
- Step 2: f(0.1, 1.1) = 1.2. y₂ = 1.1 + 0.1 × 1.2 = 1.22 at x = 0.2.
The exact value is y(0.2) = 2e^0.2 − 0.2 − 1 = 1.24281 (to 5 decimal places). The error is about 0.0228. Euler's method underestimates here because the solution curve bends upwards (it is convex), and each tangent line lies below the curve.
How Accurate Is Euler's Method?
Each step makes a small error of order h² (the local error), and over a fixed interval there are about 1/h steps, so the total (global) error is of order h. Euler's method is called a first-order method: halving the step size roughly halves the error.
For our test problem, estimating y(1) (exact value 3.43656):
- h = 0.1 (10 steps): 3.18748, error 0.249.
- h = 0.05 (20 steps): 3.30660, error 0.130.
- h = 0.025 (40 steps): 3.37013, error 0.066.
Each halving of h roughly halves the error, as predicted. To get 4 correct decimal places you would need thousands of steps, which is why better methods are used in practice.
Improved Euler (Heun's) Method
A simple improvement averages the gradient at the start and the predicted end of the step:
- Predict: p = yₙ + h f(xₙ, yₙ).
- Correct: yₙ₊₁ = yₙ + (h/2)[f(xₙ, yₙ) + f(xₙ₊₁, p)].
For the first step of our problem: p = 1.1, then y₁ = 1 + 0.05 × (1 + 1.2) = 1.11, compared with the exact 1.11034. This is a second-order method: halving h divides the error by about 4.
The Fourth-Order Runge–Kutta Method (RK4)
RK4 samples the gradient four times within each step and takes a weighted average:
- k₁ = f(xₙ, yₙ)
- k₂ = f(xₙ + h/2, yₙ + (h/2)k₁)
- k₃ = f(xₙ + h/2, yₙ + (h/2)k₂)
- k₄ = f(xₙ + h, yₙ + h k₃)
- yₙ₊₁ = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)
k₁ is the gradient at the start, k₂ and k₃ are two estimates at the midpoint, and k₄ is at the end. The weights 1, 2, 2, 1 are the same as in Simpson's rule for integration.
Example 2. Use RK4 with h = 0.1 to estimate y(0.1) for dy/dx = x + y, y(0) = 1.
- k₁ = f(0, 1) = 1
- k₂ = f(0.05, 1 + 0.05 × 1) = f(0.05, 1.05) = 1.1
- k₃ = f(0.05, 1 + 0.05 × 1.1) = f(0.05, 1.055) = 1.105
- k₄ = f(0.1, 1 + 0.1 × 1.105) = f(0.1, 1.1105) = 1.2105
- y₁ = 1 + (0.1/6)(1 + 2.2 + 2.21 + 1.2105) = 1 + (0.1/6)(6.6205) = 1.1103417
The exact value is 1.1103418. A single RK4 step is correct to 6 decimal places, while Euler's step gave 1.1.
A second step gives y(0.2) ≈ 1.2428051, against the exact 1.2428055.
Comparing Euler's Method and Runge–Kutta
RK4 is a fourth-order method: halving h divides the error by about 16. For y(1) in our test problem:
- h = 0.1: RK4 gives 3.4365595, error about 4.2 × 10⁻⁶ (Euler's error: 0.249).
- h = 0.05: error about 2.7 × 10⁻⁷.
- h = 0.025: error about 1.7 × 10⁻⁸.
RK4 does four function evaluations per step, against Euler's one, but its accuracy is so much higher that it is far cheaper for the same precision. With just 10 steps, RK4 beats Euler's method with tens of thousands of steps.
Euler's Method and RK4 in Python
The arithmetic is repetitive, so it is ideal for a computer:
import math
def f(x, y):
return x + y
def euler(f, x, y, h, steps):
for _ in range(steps):
y = y + h * f(x, y)
x = x + h
return y
def rk4(f, x, y, h, steps):
for _ in range(steps):
k1 = f(x, y)
k2 = f(x + h / 2, y + h / 2 * k1)
k3 = f(x + h / 2, y + h / 2 * k2)
k4 = f(x + h, y + h * k3)
y = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
x = x + h
return y
exact = 2 * math.exp(1) - 2
for h in [0.1, 0.05]:
n = round(1 / h)
e = euler(f, 0, 1, h, n)
r = rk4(f, 0, 1, h, n)
print(f"h={h}: Euler {e:.6f} (error {exact - e:.2e}), RK4 {r:.8f} (error {exact - r:.2e})")
Output:
h=0.1: Euler 3.187485 (error 2.49e-01), RK4 3.43655949 (error 4.17e-06) h=0.05: Euler 3.306595 (error 1.30e-01), RK4 3.43656339 (error 2.72e-07)
Systems and Higher-Order Equations
A second-order equation such as d²y/dx² = −y can be written as two first-order equations by letting v = dy/dx: then dy/dx = v and dv/dx = −y. Euler's method and RK4 apply to both equations together, updating y and v at each step. This is how computers simulate pendulums, planetary orbits and electrical circuits.
Choosing a Step Size
- Too large a step gives inaccurate or even unstable results; Euler's method applied to a rapidly decaying solution can blow up if h is too big.
- Too small a step wastes computation and lets rounding errors build up.
- Professional solvers use adaptive step sizes: they estimate the error at each step (for example by comparing two methods of different order) and shrink or grow h to keep it within a tolerance.
Summary
- Euler's method: yₙ₊₁ = yₙ + h f(xₙ, yₙ). It is simple but only first-order accurate.
- Heun's method averages two gradients and is second-order.
- RK4 combines four gradient samples with weights 1, 2, 2, 1 and is fourth-order: halving h divides the error by about 16.
- Higher-order equations are solved by rewriting them as systems of first-order equations.