Runge Kutta algorithm for numerical solution of ordinary differential equations

The solution of an ordinary differential equation (ODE) is a function. The function is of a form that it satisfies the relationship between itself and its derivatives as presented by the ODE.

Numerically solving the ODE consists of evaluating the derivative and computing the value of the function at a nearby point. This process is then repeated by evaluating the derivative at the new point and then evaluating the function at the next nearby point. The simplest finite difference algorithm for accomplishing this is called Euler's Method.

$\frac{dy}{dt}=f\left(t,y\right),y\left(t_0\right)=y_0\to \frac{y_1-y_0}{t_1-t_0}=f\left(t_0,y_0\right)$ , thus $y_1=y_0+hf\left(t_0,y_0\right),h=t_1-t_0$

Applied over and over, you get $y_{n+1}=y_n+hf\left(t_n,y_n\right),h=t_{n+1}-t_n,n=0,1,2,\dots$

While Euler's Method is appealing for its simplicity, the algorithm suffers from a low order of accuracy, it is first order. If $h$ is halved, the error in the solution is halved. This is not acceptable for real calculations.

The most common algorithm used in practice to solve ODEs is the fourth order Runge-Kutta Method because when $h$ is halved, the error is reduced by a factor of $2^4$ or 16. The fourth order Runge-Kutta Method uses a more complicated algorithm for computing the derivative of the unknown function. This leads to a very accurate result in most cases.

$\frac{dy}{dt}=f\left(t,y\right),y\left(t_0\right)=y_0,\to y_{n+1}=y_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$ where

$k_1=f\left(t_n,y_n\right),k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h{k}_1}{2}\right),k_3=f\left(t_n+\frac{h}{2},y_n+\frac{h{k}_2}{2}\right),k_4=f\left(t_{n+1},y_n+h{k}_3\right)$ and $t_n=t_0+nh$ with $h$ as the selected step size.

While this Runge-Kutta algorithm would be tedious to do on a hand-calculator, it is relatively easy to program this algorithm in many programming languages, including MATLAB. Variants of this algorithm are used in the dsolve(numerical) command in Maple and the ode45 function in MATLAB. In addition to the above steps, these computing tools estimate the error and select the value of $h$ to reduce the error. They can be used to approximate to a user specified error level value, such as ${10}^{-6}$