Runge–Kutta 4th Order (RK4) Method

Background and History of Runge–Kutta Methods

The Runge–Kutta methods are named after Carl Runge and Wilhelm (Martin) Kutta, who
developed these techniques around 1900 for solving ordinary differential equations. Early work
began with Runge (1895) and Heun (1900), but the classical 4th-order RK formula was
introduced by Kutta in 1901. In fact, Kutta’s formula is often called the ODE analog of
Simpson’s rule for numerical integration. Over the following decades (e.g. by Nyström in the
1920s) the RK family was extended, but the fourth-order scheme (RK4) remains the most widely
used ODE solver.
The Runge–Kutta methods are numerical techniques used to solve ordinary differential equations
(ODEs) of the form:
dy
=f (x , y ), y (x 0 )= y 0
dx

 Developed around 1895–1901 by Carl Runge and later extended by Martin Wilhelm
Kutta.
 Their goal was to improve upon Euler’s method, which is simple but inaccurate.
 The 4th order Runge–Kutta method (RK4) became the most popular because:
o High accuracy

o No need for higher derivatives

o Good balance between accuracy and computation cost

RK4 has a local truncation error of order O ( h5 ) and global error of order O(h 4 ).

Intuition / Idea Behind the Derivation

Why Euler’s Method is Not Enough:

Euler’s method uses only the slope at the beginning of the interval:
y n +1= y n +hf (x n , y n)

This can be inaccurate if the slope changes rapidly.

Core Idea of RK4:
Instead of using one slope, RK4 takes a weighted average of four slopes:
1. Slope at the beginning
2. Slope at the midpoint (using slope 1)
3. Slope at the midpoint (using slope 2)
4. Slope at the end
This mimics the Taylor series expansion up to 4th order without explicitly computing
derivatives.

How the RK4 Formula Is Obtained

We start with the ordinary differential equation
dy /dx = f (x, y)
with initial value y(xₙ) = yₙ.
Our goal is to find yₙ₊₁ ≈ y (xₙ + h).
1. Start from the Exact Solution (Taylor Expansion)
The exact solution can be expanded using Taylor series about xₙ:
y (xₙ + h) = yₙ
 h y′(xₙ)
 h²/2 y″(xₙ)
 h³/6 y‴(xₙ)
 h⁴/24 y⁽⁴⁾(xₙ)
 higher-order terms
Since
y′ = f (x, y),
the higher derivatives y″, y‴, y⁽⁴⁾ depend on x and y and become very complicated to compute.
Runge–Kutta avoids computing these derivatives directly.

How the Taylor Expansion of y (xₙ + h) Is Obtained:

What Is Taylor Series?
The Taylor series of a function y(x) about a point x = xₙ is:
y(x) = y(xₙ)
 (x − xₙ) y′(xₙ)
 (x − xₙ)² / 2! · y″(xₙ)
 (x − xₙ)³ / 3! · y‴(xₙ)
 (x − xₙ)⁴ / 4! · y⁽⁴⁾(xₙ)
 …
This formula comes directly from calculus and is used to approximate a function near a point.

Substitute x = xₙ + h
We want the value of y at xₙ + h, so we substitute:
x − xₙ = h
Now the Taylor series becomes:
Y (xₙ + h) = y(xₙ)
 h y′(xₙ)
 h² / 2! · y″(xₙ)
 h³ / 3! · y‴(xₙ)
 h⁴ / 4! · y⁽⁴⁾(xₙ)
 higher-order terms

Replace Factorials with Numbers

2! = 2
3! = 6
4! = 24
So the expansion becomes:
y (xₙ + h) = yₙ
 h y′(xₙ)
 h² / 2 · y″(xₙ)
  h³ / 6 · y‴(xₙ)
 h⁴ / 24 · y⁽⁴⁾(xₙ)
 higher-order terms

2. Replace Derivatives with Slope Evaluations

Instead of computing y″, y‴, etc., RK4 uses function evaluations (slopes) at different points.
First slope (start of interval)
k₁ = h f (xₙ, yₙ)
This represents the slope at xₙ.

Second slope (midpoint using k₁)

We estimate the value of y at the midpoint:
y ≈ yₙ + k₁/2
x ≈ xₙ + h/2
So,
k₂ = h f (xₙ + h/2, yₙ + k₁/2)

Third slope (refined midpoint)

Use k₂ to improve the midpoint estimate:
y ≈ yₙ + k₂/2
x ≈ xₙ + h/2
So,
k₃ = h f (xₙ + h/2, yₙ + k₂/2)

Fourth slope (end of interval)

Estimate y at the end of the interval:
y ≈ yₙ + k₃
x ≈ xₙ + h
So,
k₄ = h f (xₙ + h, yₙ + k₃)

3. Combine the Slopes

We now assume the numerical solution has the form:
yₙ₊₁ = yₙ + a k₁ + b k₂ + c k₃ + d k₄
where a, b, c, d are constants to be chosen.

4. Matching with Taylor Series

When k₁, k₂, k₃, k₄ are expanded using Taylor series and substituted into the above equation,
we compare coefficients with the exact Taylor expansion of y (xₙ + h).
By matching terms up to h⁴, we obtain:
a = 1/6
b = 1/3
c = 1/3
d = 1/6
Expand Each k-Term Using Taylor Series
We start from:
k₁ = h f(xₙ, yₙ) = h y′
So:
k₁ = h y′
Expand k₂
k₂ = h f(xₙ + h/2, yₙ + k₁/2)
Using Taylor expansion around (xₙ, yₙ):
k₂ = h [ y′ + (h/2) y″ + (h²/8) y‴ + (h³/48) y⁽⁴⁾ + … ]
Calculation of k2 expansion:
The Taylor expansion of y′(x) about x = xₙ is:
y′(xₙ + Δx)
= y′(xₙ)
 Δx · y″(xₙ)
 Δx²/2 · y‴(xₙ)
 Δx³/6 · y⁽⁴⁾(xₙ)
 …
Here, Δx = h/2.
Substitute Δx = h/2
y′(xₙ + h/2) = y′
 (h/2) y″
 (h/2)² / 2 · y‴
 (h/2)³ / 6 · y⁽⁴⁾
 …
Simplify Each Term:
(h/2)² / 2 = h² / 8
(h/2)³ / 6 = h³ / 48…

Expand k₃
k₃ = h f(xₙ + h/2, yₙ + k₂/2)
Its expansion turns out to be:
k₃ = h [ y′ + (h/2) y″ + (h²/8) y‴ + (h³/48) y⁽⁴⁾ + … ]
(For RK4 accuracy, k₂ and k₃ contribute similarly.)
Expand k₄
k₄ = h f(xₙ + h, yₙ + k₃)
Taylor expansion gives:
k₄ = h [ y′ + h y″ + h²/2 y‴ + h³/6 y⁽⁴⁾ + … ]

Substitute into the RK Formula

Substitute k₁, k₂, k₃, k₄ into:
yₙ₊₁ = yₙ + a k₁ + b k₂ + c k₃ + d k₄
Collect terms by powers of h:
Coefficient of h y′
h y′ (a + b + c + d)
This must equal h y′ from the Taylor expansion.
So:
a + b + c + d = 1 …(1)
Coefficient of h² y″
h² y″ (b/2 + c/2 + d)
This must equal h²/2 y″.
So:
b/2 + c/2 + d = 1/2 …(2)
Coefficient of h³ y‴
h³ y‴ (b/8 + c/8 + d/2)
This must equal h³/6 y‴.
So:
b/8 + c/8 + d/2 = 1/6 …(3)
Coefficient of h⁴ y⁽⁴⁾
h⁴ y⁽⁴⁾ (b/48 + c/48 + d/6)
This must equal h⁴/24 y⁽⁴⁾.
So:
b/48 + c/48 + d/6 = 1/24 …(4)
4. Solve the System of Equations
From symmetry of the method:
b=c
Now solve the equations:
From (1):
a + 2b + d = 1
From (2):
b + d = 1/2
From (3):
b/4 + d/2 = 1/6
Solving gives:
a = 1/6
b = 1/3
c = 1/3
d = 1/6

5. Final RK4 Formula

Substituting the coefficients:
1 1 1 1
y n +1= y n + k 1+ k 2+ k 3 + k 4
6 3 3 6

So:
1 2 2 1
y n +1= y n + k 1+ k 2 + k 3+ k 4
6 6 6 6
1
Factor out
6
1
y n +1= y n + ( k 1+ 2k 2+ 2k 3+ k 4 )
6

Final RK4 Formula

k 1+ 2 k 2+ 2 k 3 +k 4
y n +1= y n +
6
xₙ₊₁ = xₙ + h

Interpretation of Slopes
 k₁ → slope at the beginning of the interval
 k₂, k₃ → slopes at the midpoint
  k₄ → slope at the end of the interval
The ratio of weights is 1 : 2 : 2 : 1, which gives RK4 its high accuracy.

Examples
Problem 1:
Solve the initial value problem:
dy/dx = x + y
y(0) = 1
Use RK4 with step size h = 0.1 to find y(0.1).

Step 1: Initial values

x₀ = 0
y₀ = 1
h = 0.1

Step 2: Compute k₁
k₁ = 0.1 × (0 + 1)
k₁ = 0.1

Step 3: Compute k₂
x₀ + h/2 = 0.05
y₀ + k₁/2 = 1 + 0.05 = 1.05
k₂ = 0.1 × (0.05 + 1.05)
k₂ = 0.11

Step 4: Compute k₃
x₀ + h/2 = 0.05
y₀ + k₂/2 = 1 + 0.055 = 1.055
k₃ = 0.1 × (0.05 + 1.055)
k₃ = 0.1105

Step 5: Compute k₄
x₀ + h = 0.1
y₀ + k₃ = 1 + 0.1105 = 1.1105
k₄ = 0.1 × (0.1 + 1.1105)
k₄ = 0.12105

Step 6: Compute y₁
y₁ = 1 + (0.1 + 2(0.11) + 2(0.1105) + 0.12105) / 6
y₁ = 1 + 0.66205 / 6
y₁ ≈ 1.11034

Final Answer
y(0.1) ≈ 1.11034
(The exact solution is approximately 1.1103418, showing that RK4 gives a very accurate result.)

Problem 2:
Solve the initial value problem:
dy/dx = x − y
y(0) = 1
Use RK4 with step size h = 0.1 to find y(0.1).

Step 1: Initial values

x₀ = 0
y₀ = 1
h = 0.1
Step 2: Compute k₁
k₁ = 0.1 × (0 − 1)
k₁ = −0.1
Step 3: Compute k₂
x₀ + h/2 = 0.05
y₀ + k₁/2 = 1 − 0.05 = 0.95
k₂ = 0.1 × (0.05 − 0.95)
k₂ = −0.09
Step 4: Compute k₃
x₀ + h/2 = 0.05
y₀ + k₂/2 = 1 − 0.045 = 0.955
k₃ = 0.1 × (0.05 − 0.955)
k₃ = −0.0905
Step 5: Compute k₄
x₀ + h = 0.1
y₀ + k₃ = 1 − 0.0905 = 0.9095
k₄ = 0.1 × (0.1 − 0.9095)
k₄ = −0.08095
Step 6: Compute y₁
y₁ = 1 + (−0.1 + 2(−0.09) + 2(−0.0905) + (−0.08095)) / 6
y₁ = 1 − 0.542 / 6
y₁ ≈ 0.90967
Final Answer
y(0.1) ≈ 0.90967

Problem 3:
Solve the initial value problem:
dy/dx = x² + y
y(0) = 1
Use RK4 with step size h = 0.1 to find y(0.1).

Step 1: Initial values

x₀ = 0
y₀ = 1
h = 0.1
Step 2: Compute k₁
k₁ = 0.1 × (0² + 1)
k₁ = 0.1
Step 3: Compute k₂
x₀ + h/2 = 0.05
y₀ + k₁/2 = 1 + 0.05 = 1.05
k₂ = 0.1 × (0.05² + 1.05)
k₂ = 0.10525
Step 4: Compute k₃
x₀ + h/2 = 0.05
y₀ + k₂/2 = 1 + 0.052625 = 1.052625
k₃ = 0.1 × (0.05² + 1.052625)
k₃ = 0.10551
Step 5: Compute k₄
x₀ + h = 0.1
y₀ + k₃ = 1 + 0.10551 = 1.10551
k₄ = 0.1 × (0.1² + 1.10551)
k₄ = 0.111551
Step 6: Compute y₁
y₁ = 1 + (0.1 + 2(0.10525) + 2(0.10551) + 0.111551) / 6
y₁ = 1 + 0.632071 / 6
y₁ ≈ 1.10535
Final Answer
y(0.1) ≈ 1.10535

Problem 4 :
Solve the initial value problem:
dy/dx = y − x²
y(0) = 1
Use RK4 with step size h = 0.1 to find y(0.1).
Step 1: Initial values
x₀ = 0
y₀ = 1
h = 0.1
Step 2: Compute k₁
k₁ = 0.1 × (1 − 0²)
k₁ = 0.1
Step 3: Compute k₂
x₀ + h/2 = 0.05
y₀ + k₁/2 = 1 + 0.05 = 1.05
k₂ = 0.1 × (1.05 − 0.05²)
k₂ = 0.10475
Step 4: Compute k₃
x₀ + h/2 = 0.05
y₀ + k₂/2 = 1 + 0.052375 = 1.052375
k₃ = 0.1 × (1.052375 − 0.05²)
k₃ = 0.10499
Step 5: Compute k₄
x₀ + h = 0.1
y₀ + k₃ = 1 + 0.10499 = 1.10499
k₄ = 0.1 × (1.10499 − 0.1²)
k₄ = 0.109499
Step 6: Compute y₁
y₁ = 1 + (0.1 + 2(0.10475) + 2(0.10499) + 0.109499) / 6
y₁ = 1 + 0.628979 / 6
y₁ ≈ 1.10483
Final Answer
y(0.1) ≈ 1.10483