Runge-Kutta Methods

In the forward Euler method, we used the information on the slope or the derivative of y at the given time step to extrapolate the solution to the next time-step. The LTE for the method is O(h²), resulting in a first order numerical technique. Runge-Kutta methods are a class of methods which judiciously uses the information on the 'slope' at more than one point to extrapolate the solution to the future time step. Let's discuss first the derivation of the second order RK method where the LTE is O(h³).

Given the IVP of Eq. 6, and a time step h, and the solution y_n at the nth time step, let's say that we wish to compute y_n+1 in the following fashion:

$$k_2 = hf(y_n+\beta k_1, t_n + \alpha h)$$

y_n+1 = y_n + ak₁ + bk₂, (12)

where the constants $\alpha$, $\beta$, a and b have to be evaluated so that the resulting method has a LTE O(h³). Note that if k₂=0 and a=1, then Eq. 13 reduces to the forward Euler method.

Now, let's write down the Taylor series expansion of y in the neighborhood of t_n correct to the h² term i.e.,

$$y(t_{n+1}) = y(t_n) + h \frac{dy}{dt}\vert _{t_n} + \frac {h^2}{2} \frac{d^2y}{dt^2}\vert _{t_n} + {\mbox{O}}(h^3).$$ (13)

However, we know from the IVP (Eq. 6) that dy/dt = f(y,t) so that

$$\frac{d^2y}{dt^2} = \frac{df(y,t)}{dt} = \frac{\partial f}{\p... ...rac{\partial f}{\partial t} + f \frac{\partial f}{\partial y}.$$ (14)

So from the above analysis, i.e., Eqs. 14 and 15, we get

$$y_{n+1} = y_n + h f(y_n,t_n) + \frac{h^2}{2} \left[ \frac{\pa... ...ac{\partial f}{\partial y}\right](y_n,t_n) + {\mbox{O}}(h^3).$$ (15)

However, the term k₂ in the proposed RK method of Eq. 13 can be expanded correct to O(h³) as

$$k_2 = hf(y_n+\beta k_1, t_n + \alpha h)$$

$$= h\left( f(y_n,t_n) + \alpha h\frac{\partial f}{\partial t} (y_n... ... + \beta k_1 \frac{\partial f}{\partial t} (y_n, t_n)\right) + {\mbox{O}}(h^3).$$ (16)

Now, substituting for k₂ from Eq. 17 in Eq. 13, we get

$$y_{n+1} = y_n + (a+b) h f(y_n,t_n) + bh^2(\alpha \frac {\part... ...f \frac {\partial f}{\partial y})(y_n,t_n) + {\mbox{O}}(h^3).$$ (17)

Comparing the terms with identical coefficients in Eqs. 16 and 18 gives us the following system of equations to determine the constants:

$$\alpha b = \frac{1}{2}$$

$$\beta b = \frac{1}{2}.$$ (18)

There are infinitely many choices of a, b, $\alpha$ and $\beta$ which satisfy Eq. 19, we can choose for instance $\alpha = \beta = 1$ and a=b=1/2. With this choice, we have the classical second order accurate Runge-Kutta method (RK2) which is summarized as follows.

$$y_{n+1} = y_n + (k_1 + k_2)/2, \:\:\: {\mbox{Second Order Runge-Kutta Method (RK2).}}$$ (19)

In a similar fashion Runge-Kutta methods of higher order can be developed. One of the most widely used methods for the solution of IVPs is the fourth order Runge-Kutta (RK4) technique. The LTE of this method is order h⁵. The method is given below.

$$k_3 = hf(y_n+k_2/2, t_n + h/2) \:\: \: {\mbox{Fourth Order Runge-Kutta Method (RK4).}}$$ k₄ = h(y_n+k₃, t_n + h)

y_n+1 = y_n + (k₁ + 2k₂ + 2k₃ + k₄)/6. (20)

Note that the RK methods are explicit techniques, hence they are only conditionally stable.