Runge–Kutta–Fehlberg method

Any Runge–Kutta method is uniquely identified by its Butcher tableau. The embedded pair proposed by Fehlberg[2]

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.

This shows the computational time in real time used during a 3-body simulation evolved with the Runge-Kutta-Fehlberg method. Most of the computer time is spent when the bodies pass close by and are susceptible to numerical error.

The coefficients found by Fehlberg for Formula 1 (derivation with his parameter α2=1/3) are given in the able below, using array indexing of base 1 instead of base 0 to be compatible with most computer languages:

Fehlberg[2] outlines a solution to solving a system of n differential equations of the form:

${\displaystyle {\operatorname {d} \!y_{i} \over \operatorname {d} \!x}=f_{i}(x,y_{1},y_{2},...,y_{n}),i=1,2,...,y_{n}}$

to iterative solve for

${\displaystyle y_{i}(x+h),i=1,2,...,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, ${\textstyle h}$, and an estimated slope specified by function f on the right-hand side of the differential equation.

${\displaystyle k_{1}=h\cdot f(x+A(1)\cdot h,y)}$

${\displaystyle k_{2}=h\cdot f(x+A(2)\cdot h,y+B(2,1)\cdot k_{1})}$

${\displaystyle k_{3}=h\cdot f(x+A(3)\cdot h,y+B(3,1)\cdot k_{1}+B(3,2)\cdot k_{2})}$

${\displaystyle 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})}$

${\displaystyle 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})}$

${\displaystyle 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})}$

Then the weighted average is:

${\displaystyle 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:

${\displaystyle TE=|CT(1)\cdot k_{1}+CT(2)\cdot k_{2}+CT(3)\cdot k_{3}+CT(4)\cdot k_{4}+CT(5)\cdot k_{5}+CT(6)\cdot k_{6}|}$

At the completion of the step, a new stepsize is calculated:

${\displaystyle h_{new}=0.9\cdot h\cdot \left({\frac {\epsilon }{|TE|}}\right)^{1/5}}$

If ${\textstyle |TE|>\epsilon }$, then replace ${\textstyle h}$ with ${\textstyle h_{new}}$ and repeat the step. If ${\textstyle |TE|\leqslant \epsilon }$, then the step is completed. Replace ${\textstyle h}$ with ${\textstyle h_{new}}$ for the next step.

The coefficients found by Fehlberg for Formula 2 (derivation with his parameter α2=3/8) are given in the able below, using array indexing of base 1 instead of base 0 to be compatible with most computer languages:

In another table in Fehlberg,[2] coefficients for an RKF4(5) derived by D. Sarafyan are given:

• List of Runge–Kutta methods
• Numerical methods for ordinary differential equations
• Runge–Kutta methods

• Erwin Fehlberg (1969). Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems . NASA Technical Report 315. https://ntrs.nasa.gov/api/citations/19690021375/downloads/19690021375.pdf
• Erwin Fehlberg (1968) Classical fifth-, sixth-, seventh-, and eighth-order Runge-Jutta formulas with stepsize control. NASA Technical Report 287. https://ntrs.nasa.gov/api/citations/19680027281/downloads/19680027281.pdf
• Erwin Fehlberg (1970) Some experimental results concerning the error propagation in Runge-Kutta type integration formulas. NASA Technical Report R-352. https://ntrs.nasa.gov/api/citations/19700031412/downloads/19700031412.pdf
• Erwin Fehlberg (1970). "Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Wärmeleitungsprobleme," Computing (Arch. Elektron. Rechnen), vol. 6, pp. 61–71. doi:10.1007/BF02241732
• Ernst Hairer, Syvert Nørsett, and Gerhard Wanner (1993). Solving Ordinary Differential Equations I: Nonstiff Problems, second edition, Springer-Verlag, Berlin. ISBN 3-540-56670-8.
• Diran Sarafyan (1966) Error Estimation for Runge-Kutta Methods Through Pseudo-Iterative Formulas. Technical Report No. 14, Louisiana State University in New Orleans, May 1966.

• Simos, T. E. (1993). A Runge-Kutta Fehlberg method with phase-lag of order infinity for initial-value problems with oscillating solution. Computers & Mathematics with Applications, 25(6), 95-101.
• Handapangoda, C. C., Premaratne, M., Yeo, L., & Friend, J. (2008). Laguerre Runge-Kutta-Fehlberg Method for Simulating Laser Pulse Propagation in Biological Tissue. IEEE Journal of Selected Topics in Quantum Electronics, 14(1), 105-112.
• Paul, S., Mondal, S. P., & Bhattacharya, P. (2016). Numerical solution of Lotka Volterra prey predator model by using Runge–Kutta–Fehlberg method and Laplace Adomian decomposition method. Alexandria Engineering Journal, 55(1), 613-617.
• Filiz, A. (2014). Numerical solution of linear Volterra integro-differential equation using Runge–Kutta–Fehlberg method. Applied and Computational Mathematics, 3(1), 9-14.
• Simos, T. E. (1995). Modified Runge–Kutta–Fehlberg methods for periodic initial-value problems. Japan journal of industrial and applied mathematics, 12(1), 109.
• Sarafyan, D. (1994) Approximate Solution of Ordinary Differential Equations and Their Systems Through Discrete and Continuous Embedded Runge-Kutta Formulae and Upgrading Their Order, Computers Math. Applic. Vol. 28, No. 10-12, pp. 353–384, 1994 https://core.ac.uk/download/pdf/82540775.pdf