Overview
By extending the Taylor expansion of the function f up to order 5, a reasoning similar to that used for second-order Runge–Kutta methods leads to a nonlinear system of eight equations containing ten unknowns (see Scheid, ref. [32]). The final result is the classical fourth-order Runge–Kutta method, which is a highly useful and widely applied tool.
Runge–Kutta Order 4 Algorithm
Given: step size h, initial condition (t₀, y₀), and a maximum number of iterations N.
For 0 ≤ n ≤ N:
k₁ = h f(tₙ, yₙ)
k₂ = h f(tₙ + h/2 , yₙ + k₁/2)
k₃ = h f(tₙ + h/2 , yₙ + k₂/2)
k₄ = h f(tₙ + h , yₙ + k₃)
yₙ₊₁ = yₙ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄)
tₙ₊₁ = tₙ + h
Write tₙ₊₁ and yₙ₊₁.
Remark
The Runge–Kutta method of order 4 is very frequently used because of its high accuracy, which is clearly illustrated in the following example.
Example
Consider again the differential equation:
y′(t) = −y(t) + t + 1 , y(0) = 1
We evaluate the constants kᵢ. For the first iteration (h = 0.1):
k₁ = 0.1 f(0, 1) = 0.1(−1 + 1) = 0
k₂ = 0.1 f(0.05, 1) = 0.1(−1 + 1.05) = 0.005
k₃ = 0.1 f(0.05, 1.0025) = 0.1(−1.0025 + 1.05) = 0.00475
k₄ = 0.1 f(0.1, 1.00475) = 0.1(−1.00475 + 1.1) = 0.009525
y₁ = 1 + (1/6)(0 + 2·0.005 + 2·0.00475 + 0.009525)
= 1.0048375
Second iteration:
k₁ = 0.1 f(0.1, 1.0048375) = 0.00951625
k₂ = 0.1 f(0.15, 1.009595625) = 0.014040438
k₃ = 0.1 f(0.15, 1.011857719) = 0.0138142281
k₄ = 0.1 f(0.2, 1.018651728) = 0.0181348272
y₂ = 1.0048375 + (1/6)(k₁ + 2k₂ + 2k₃ + k₄)
= 1.0187309014
The following table compares the numerical and exact solutions, and gives the absolute error.
| tᵢ | y(tᵢ) | yᵢ | |y(tᵢ) − yᵢ| |
|---|---|---|---|
| 0.0 | 1.0 | 1.0 | 0.0 |
| 0.1 | 1.0048374180 | 1.0048375000 | 0.819×10⁻⁷ |
| 0.2 | 1.0187307798 | 1.0187309014 | 0.148×10⁻⁶ |
| 0.3 | 1.0408182207 | 1.0408184220 | 0.210×10⁻⁶ |
| 0.4 | 1.0703200460 | 1.0703202889 | 0.242×10⁻⁶ |
| 0.5 | 1.1065306597 | 1.1065309344 | 0.274×10⁻⁶ |
| 0.6 | 1.1488116361 | 1.1488119343 | 0.298×10⁻⁶ |
| 0.7 | 1.1965853034 | 1.1965856186 | 0.314×10⁻⁶ |
| 0.8 | 1.2493289641 | 1.2493292897 | 0.325×10⁻⁶ |
| 0.9 | 1.3065696598 | 1.3065799912 | 0.331×10⁻⁶ |
| 1.0 | 1.3678794412 | 1.3678797744 | 0.333×10⁻⁶ |
The error is on the order of 10⁻⁶, which compares very favorably with lower-order methods. A slight increase in error over successive iterations is also observed, showing once again the propagation of error from one step to the next.