# Runge–Kutta–Fehlberg method

In mathematics, the Runge–Kutta–Fehlberg method (or Fehlberg method) is an algorithm in numerical analysis for the numerical solution of ordinary differential equations. It was developed by the German mathematician Erwin Fehlberg and is based on the large class of Runge–Kutta methods.
The novelty of Fehlberg's method is that it is an embedded method from the Runge–Kutta family, meaning that it reuses the same intermediate calculations to produce two estimates of different accuracy, allowing for automatic error estimation. The method presented in Fehlberg's 1969 paper has been dubbed the RKF45 method, and is a method of order O(h$^4$) with an error estimator of order O(h$^5$). By performing one extra calculation, the error in the solution can be estimated and controlled by using the higher-order embedded method that allows for an adaptive stepsize to be determined automatically.

## Butcher tableau for Fehlberg's 4(5) method

Any Runge–Kutta method is uniquely identified by its Butcher tableau. The embedded pair proposed by Fehlberg
| 0 |
| 1/4 | 1/4 |
| 3/8 | 3/32 | 9/32 |
| 12/13 | 1932/2197 | −7200/2197 | 7296/2197 |
| 1 | 439/216 | −8 | 3680/513 | −845/4104 |
| 1/2 | −8/27 | 2 | −3544/2565 | 1859/4104 | −11/40 | |
| | 16/135 | 0 | 6656/12825 | 28561/56430 | −9/50 | 2/55 |
| | 25/216 | 0 | 1408/2565 | 2197/4104 | −1/5 | 0 |
The first row of coefficients at the bottom of the table gives the fifth-order accurate method, and the second row gives the fourth-order accurate method.

## Implementing an RK4(5) Algorithm

The coefficients found by Fehlberg for Formula 1 (derivation with his parameter α$_2$=1/3) are given in the table below, using array indexing of base 1 instead of base 0 to be compatible with most computer languages:
COEFFICIENTS FOR RK4(5), FORMULA 1 Table II in Fehlberg :
| K | A(K) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | C(K) | CH(K) | CT(K) |
| K | A(K) | L=1 | L=2 | L=3 | L=4 | L=5 | C(K) | CH(K) | CT(K) |
| 1 | 0 | | | | | | 1/9 | 47/450 | 1/150 |
| 2 | 2/9 | 2/9 | | | | | 0 | 0 | 0 |
| 3 | 1/3 | 1/12 | 1/4 | | | | 9/20 | 12/25 | -3/100 |
| 4 | 3/4 | 69/128 | -243/128 | 135/64 | | | 16/45 | 32/225 | 16/75 |
| 5 | 1 | -17/12 | 27/4 | -27/5 | 16/15 | | 1/12 | 1/30 | 1/20 |
| 6 | 5/6 | 65/432 | -5/16 | 13/16 | 4/27 | 5/144 | | 6/25 | -6/25 |
Fehlberg outlines a solution to solving a system of n differential equations of the form: ${\frac {dy_{i}}{dx}}=f_{i}(x,y_{1},y_{2},\ldots ,y_{n}),i=1,2,\ldots ,n$ to iterative solve for $y_{i}(x+h),i=1,2,\ldots ,n$ where h is an adaptive stepsize to be determined algorithmically:
The solution is the weighted average of six increments, where each increment is the product of the size of the interval, h{\textstyle h}, and an estimated slope specified by function f on the right-hand side of the differential equation.
${\begin{aligned}k_{1}&=h\cdot f(x+A(1)\cdot h,y)\\k_{2}&=h\cdot f(x+A(2)\cdot h,y+B(2,1)\cdot k_{1})\\k_{3}&=h\cdot f(x+A(3)\cdot h,y+B(3,1)\cdot k_{1}+B(3,2)\cdot k_{2})\\k_{4}&=h\cdot f(x+A(4)\cdot h,y+B(4,1)\cdot k_{1}+B(4,2)\cdot k_{2}+B(4,3)\cdot k_{3})\\k_{5}&=h\cdot f(x+A(5)\cdot h,y+B(5,1)\cdot k_{1}+B(5,2)\cdot k_{2}+B(5,3)\cdot k_{3}+B(5,4)\cdot k_{4})\\k_{6}&=h\cdot f(x+A(6)\cdot h,y+B(6,1)\cdot k_{1}+B(6,2)\cdot k_{2}+B(6,3)\cdot k_{3}+B(6,4)\cdot k_{4}+B(6,5)\cdot k_{5})\end{aligned}}$
Then the weighted average is:
$y(x+h)=y(x)+CH(1)\cdot k_{1}+CH(2)\cdot k_{2}+CH(3)\cdot k_{3}+CH(4)\cdot k_{4}+CH(5)\cdot k_{5}+CH(6)\cdot k_{6}$
The estimate of the truncation error is: $\mathrm {TE} =\left|\mathrm {CT} (1)\cdot k_{1}+\mathrm {CT} (2)\cdot k_{2}+\mathrm {CT} (3)\cdot k_{3}+\mathrm {CT} (4)\cdot k_{4}+\mathrm {CT} (5)\cdot k_{5}+\mathrm {CT} (6)\cdot k_{6}\right|$
At the completion of the step, a new stepsize is calculated:
$h_{\text{new}}=0.9\cdot h\cdot \left({\frac {\varepsilon }{TE}}\right)^{1/5}$
If TE>ε{\textstyle \mathrm {TE} >\varepsilon }, then replace h{\textstyle h} with h$_{new}${\textstyle h_{\text{new}}} and repeat the step. If TE⩽ε{\textstyle TE\leqslant \varepsilon }, then the step is completed. Replace h{\textstyle h} with h$_{new}${\textstyle h_{\text{new}}} for the next step.
The coefficients found by Fehlberg for Formula 2 (derivation with his parameter α$_2$ = 3/8) are given in the table below, using array indexing of base 1 instead of base 0 to be compatible with most computer languages:
COEFFICIENTS FOR RK4(5), FORMULA 2 Table III in Fehlberg :
| K | A(K) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | C(K) | CH(K) | CT(K) |
| K | A(K) | L=1 | L=2 | L=3 | L=4 | L=5 | C(K) | CH(K) | CT(K) |
| 1 | 0 | | | | | | 25/216 | 16/135 | -1/360 |
| 2 | 1/4 | 1/4 | | | | | 0 | 0 | 0 |
| 3 | 3/8 | 3/32 | 9/32 | | | | 1408/2565 | 6656/12825 | 128/4275 |
| 4 | 12/13 | 1932/2197 | -7200/2197 | 7296/2197 | | | 2197/4104 | 28561/56430 | 2197/75240 |
| 5 | 1 | 439/216 | -8 | 3680/513 | -845/4104 | | -1/5 | -9/50 | -1/50 |
| 6 | 1/2 | -8/27 | 2 | -3544/2565 | 1859/4104 | -11/40 | | 2/55 | -2/55 |
In another table in Fehlberg, coefficients for an RKF4(5) derived by D. Sarafyan are given:
COEFFICIENTS FOR Sarafyan's RK4(5), Table IV in Fehlberg :
| K | A(K) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | B(K,L) | C(K) | CH(K) | CT(K) |
| K | A(K) | L=1 | L=2 | L=3 | L=4 | L=5 | C(K) | CH(K) | CT(K) |
| 1 | 0 | 0 | | | | | 1/6 | 1/24 | 1/8 |
| 2 | 1/2 | 1/2 | | | | | 0 | 0 | 0 |
| 3 | 1/2 | 1/4 | 1/4 | | | | 2/3 | 0 | 2/3 |
| 4 | 1 | 0 | -1 | 2 | | | 1/6 | 5/48 | 1/16 |
| 5 | 2/3 | 7/27 | 10/27 | 0 | 1/27 | | | 27/56 | -27/56 |
| 6 | 1/5 | 28/625 | -1/5 | 546/625 | 54/625 | -378/625 | | 125/336 | -125/336 |