Higher Order Integrator

Runge-Kutta $2^{th}$ order - Mid-Point method

The Euler and Mid-Point methods of order $O(h^2)$ and $O(h^3)$ respectively. This can be seen via a Taylor series expansion of the solution of the differential equation. Let the equation be

$$\frac{ dy}{d t } = f( y, t ).$$

The unknown function value at $t_0+h$ is obtained via an expansion of the solution about $t=t_0$, yielding

$$y(t_0+h) = y(t_0) + h \frac{ d y}{d t }(t_0) + \frac{h^2}{2} \frac{ d y^2}{d t^2 }(t_0) +O(h^3).$$

In this equation $h=dt$ is the time step. Keeping the first wo term yields the Euler method with neglected terms of order $O(h^2)$.

The Mid-Point method takes into account of the next term in the series. This requires knowing $\frac{ d y^2}{d t^2 }(t_0)$. For simplicity in the derivation, let us assume that $f(y,t)$ depends only on $y(t)$. Then

$$\frac{ d y^2}{d t^2 } = \frac{ d f}{d t } = \frac{ \parti... ...ac{ d y}{d t } = \frac{ \partial f}{\partial y } f(y) = f' f.$$

Calculating $f'$ could be complicated and expensive. Instead this term can be approximated in terms of a Taylor series expansion of $f(y)$ and substituted in the original expansion yielding error of order $O(h^3)$. Use

$$f(y_0 + \Delta y) = f(y_0) + \Delta y f'(y_0) + O(\Delta y^2)$$

Choose $\Delta y = \frac{h}{2} f(y_0)$ to yield

$$f(y_0 + \frac{h}{2} f(y_0) ) = f(y_0) + \frac{h}{2} f(y_0) f'(y_0) + O(\Delta y^2)$$

Use $f(y_0) f'(y_0) = y^{''}_0$ and multiply by $h$ to obtain

$$h f( y_0 + \frac{h}{2} f(y_0) ) - h f(y_0) = \frac{h^2}{2} y^{''}_0 + O(h^3)$$

Substituting back in the first Taylor series yields

$$y(t_0+h) = y(t_0) + h f(y_0 + \frac{h}{2} f(y_0) )$$

This is the Mid-Point method. It is also called Runge-Kutta of second order. The neglected terms are of order $O(h^3)$. The order of an integrator is one order less than the neglected term in the Taylor series.

The Runge-Kutta of second order yields the four steps we previously deduced from an analogy with the symmetric form of the finite difference formula for the first derivative of a function on a grid:

1. Calculate $f(y_0) = f(y(t_0),t_0)$
2. Calculate $y_{half} = y_0 + \frac{h}{2} f(y_0)$
3. Calculate $f(y_{half}) = f(y_0 + \frac{h}{2} f(y_0), t_0 + \frac{h}{2} )$
4. Calculate $y(t_0+h) = y(t_0) + h f(y_{half})$

Predictor-Corrector Improved Euler revisited

Predictor-Corrector Improved Euler is of the same order as the Mid-Point method. Its derivation in terms of a Taylor series expansion illustrates this fact. It also shows that different methods result from different choice in canceling higher order terms.

Runge-Kutta $4^{th}$ order

Higher order terms in the Taylor series expansion can be included in the method of solution. The most popular extension of the Runge-Kutta of order 2 is the one of order 4. This method of solution only neglects terms of order $O(h^5)$ in each time step. The formula in this case corresponds to the following steps:

1. $k_1 = h f(y_0, t_0)$

2. $k_2 = h f( y_0+k_1/2, t_0+h/2)$

3. $k_3 = h f( y_0+k_2/2, t_0+h/2)$

4. $k_4 = h f( y_0+k_3, t_0+h)$

5. $y(t_0+h) = y(t_0) + ( k_1/2 + k_2 + k_3 + k_4/2 ) / 3$

The coefficients in this formula are chosen to guaranty that the first dropped term is of order $O(h^5)$. Its use requires 4 evaluations of the RHS function $f(y(t),t)$. But the higher order of the Runge-Kutta $4^{th}$ order integrator may allow larger time steps to be used which will compensate for this higher number of RHS function calls.

The Runge-Kutta $4^{th}$ order time step approach is illustrated in the following figure.

The C implementation of a Runge-Kutta $4^{th}$ order time step is a generalization of our previous time step general functions (e.g., Euler and Mid-Point time step functions): rk4_step.c

The C implementation of a Runge-Kutta $4^{th}$ order time step with a rhs function pointer is given in: rk4_step_pointer.c

Numerical Recipes says: For many scientific users, fourth-order Runge-Kutta is not just the first word on ODE integrators, but the last word as well.

Exercise: Solve the Harmonic Oscillator using the Runge-Kutta $4^{th}$ order integrator. Vary the time step and compare the results to the Euler and Mid-Point solutions.