In the last section it was shown that using two estimates of the slope (i.e., Second Order Runge Kutta; using slopes at the beginning and midpoint of the time step, or using the slopes at the beginninng and end of the time step) gave an approximation with greater accuracy than using just a single slope (i.e., First Order Runge Kutta; using only the slope at the beginning of the interval). It seems reasonable, then, to assume that using even more estimates of the slope would result in even more accuracy. It turns out this is true, and leads to higher-order Runge-Kutta methods. Third order methods can be developed (but are not discussed here). Instead we will restrict ourselves to the much more commonly used Fourth Order Runge-Kutta technique, which uses four approximations to the slope. It is important to understand these lower order methods before starting on the fourthe order method.

If you are interested in the details of the derivation of the Fourth Order Runge-Kutta Methods, check a Numerical Methods Textbook (like *Applied Numerical Methods*, by Carnahan, Luther and Wilkes)

To review the problem at hand: we wisth to approximate the solution to a first order differential equation given by

$${{dy(t)} \over {dt}} = y'(t) = f(y(t),t), \quad \quad {\rm{with\;}} y(t_0)=y_0$$(starting from some known initial condition, *y(t₀)=y₀*). The development of the Fourth Order Runge-Kutta method closely follows those for the Second Order, and will not be covered in detail here. As with the second order technique there are many variations of the fourth order method, and they all use four approximations to the slope. We will use the following slope approximations to estimate the slope at some time *t₀* (assuming we only have an approximation to* y(t₀) *(which we call *y*(t₀)*).

& {k_1} = f({y^*}({t_0}),{t_0}) \cr

& {k_2} = f\left( {{y^*}({t_0}) + {k_1}{h \over 2},{t_0} + {h \over 2}} \right) \cr

& {k_3} = f\left( {{y^*}({t_0}) + {k_2}{h \over 2},{t_0} + {h \over 2}} \right) \cr

& {k_4} = f\left( {{y^*}({t_0}) + {k_3}h,{t_0} + h} \right) \cr} $$

Each of these slope estimates can be described verbally.

*k*is the slope at the beginning of the time step (this is the same as_{1}*k*in the first and second order methods)._{1}- If we use the slope
*k*to step halfway through the time step, then_{1}*k*is an estimate of the slope at the midpoint. This is the same as the slope,_{2}*k*, from the second order midpoint method. This slope proved to be more accurate than_{2}*k*for making new approximations for y(t)._{1} - If we use the slope
*k*to step halfway through the time step, then k_{2}_{3}is another estimate of the slope at the midpoint. - Finally, we use the slope,
*k*, to step all the way across the time step (to_{3}*t₀+h*), and*k*is an estimate of the slope at the endpoint._{4}

We then use a weighted sum of these slopes to get our final estimate of *y*(t₀+h)*

{y^*}({t_0} + h) &= {y^*}({t_0}) + {{{k_1} + 2{k_2} + 2{k_3} + {k_4}} \over 6}h = {y^*}({t_0}) + \left( {{1 \over 6}{k_1} + {1 \over 3}{k_2} + {1 \over 3}{k_3} + {1 \over 6}{k_4}} \right)h \cr

&= {y^*}({t_0}) + mh\quad \quad {\rm{where\;}}m{\rm{\;is\;a\;weighted\;average\; slope\; approximation}} \cr} $$

Note: these weights are chosen to get the desired accuracy of the approximation, just as with the second order method. The details are not given here.

The conceptual idea is similar to the second order endpoint method which used an equal weighting of the slopes at the beginning and end of the interval. Here the weighting of the midpoint slopes (*k _{2}* and

As an example consider the differential equation

$${{dy(t)} \over {dt}} + 2y(t) = 0\quad {\rm{or }} \quad {{dy(t)} \over {dt}} = - 2y(t)$$with the initial condition set as *y(0)=3*. To get from the initial value at *t=0* to an estimate at *t=h*, follow the procedure outlined below

y'(0) &= - 2y(0)\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = 0 \cr

{k_1} &= - 2y(0)\quad &{\rm{derivative}}\;{\rm{at}}\;t = 0 \cr

{y_1}\left( {{h \over 2}} \right) &= {y^*}(0) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{1}}}{\rm{)}} \cr

{k_2} &= - 2{y_1}\left( {{h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = h/2 \cr

{y_2}\left( {{h \over 2}} \right) &= {y^*}(0) + {k_2}{h \over 2}\quad &{\rm{another intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{2}}}{\rm{)}} \cr

{k_3} &= - 2{y_2}\left( {{h \over 2}} \right)\quad &{\rm{another}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = h/2 \cr

{y_3}\left( h \right) &= {y^*}(0) + {k_3}h\quad &{\rm{an}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = h\;{\rm{(using\; }}{{\rm{k}}_{\rm{3}}}{\rm{)}} \cr

{k_4} &= - 2{y_3}\left( h \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = h \cr \cr

{y^*}(h) &= y(0) + {{{k_1} + 2{k_2} + 2{k_3} + {k_4}} \over 6}h\quad &{\rm{estimate}}\;{\rm{of}}\;y(h) \cr} $$

in general to go from an estimate at *t₀ *to an estimate at *t₀+h*,

{k_1} &= - 2{y^*}({t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{1}}}{\rm{)}} \cr

{k_2} &= - 2{y_1}\left( {{h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y_2}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_2}{h \over 2}\quad &{\rm{another}}\;{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{2}}}{\rm{)}} \cr

{k_3} &= - 2{y_2}\left( {{h \over 2}} \right)\quad &{\rm{another}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y_3}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_3}h\quad &{\rm{an}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h\;{\rm{(using}}\;{{\rm{k}}_{\rm{3}}}{\rm{)}} \cr

{k_4} &= - 2{y_3}\left( {{t_0} + h} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h \cr \cr

{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {{{k_1} + 2{k_2} + 2{k_3} + {k_4}} \over 6}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$

The Fourth Order Runge-Kutta method is fairly complicated. This section of the text is an attempt to help to visualize the process; you should feel free to skip it if it already makes sense to you and go on to the example that follows.

We will use the same problem as before. We will call the initial time *t₀* and set *t₀=0*. We will then estimate a solution of the differential equation

$${{dy(t)} \over {dt}} = y'(t) = -2y(t), \quad \quad {\rm{with\;}} y(0)=3$$

after one time step, i.e., at *t=t₀+h* or *t*=h, since *t₀=0*.

Given y(0), we find the slope, *y'(0)* using the equation stated above. We call this* slope k _{1}*.

We use this slope to estimate the value of the function midway through the intervale, i.e., *y(½h)*. We call this estimate *y _{1}*,

Given y_{1}, we can use it to estimate the slope at the midpoint of the interval, *t=½h*. We call this slope *k2*.

Note that this is an estimate of the slope at t=½h and we use it to find another estimate of *y(½h)*, that we call* y _{2}*, with

Given *y _{2}*, we can use it to find another estimate the slope at

We use this estimate of the slope at *t=½h* to find an estimate of the function at the end of the interval, *y(h)*. We call this estimate *y _{3}*, with

We use *y _{3}* to find an estimate of the slope at the end of the interval,

We use all of our estimates of the slope in a weighted average that we will use to generate our final estimate for *y(h)*. We give the midpoint slope estimates twice as much weight as the endpoint estimates. We define this weighted estimate as

and use it to generate our final estimate

$${y^*}(h) = y(0) + \left( {{1 \over 6}{k_1} + {1 \over 3}{k_2} + {1 \over 3}{k_3} + {1 \over 6}{k_4}} \right)h = y(0) + m \cdot h$$The value of this final estimate for the given example is *y*(h)=2.0112*. This is quite close to the exact solution* y(h)=3e ^{-2(0.2)}=2.0110*. Note: As stated previously, we generally won't know the exact solution as we do in this case.

We can repeat this process using this estimate as our starting point to generate a solution over time.

We can use MATLAB to perform the calculation described above.A simple loop accomplishes this:

% Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta y0 = 3; % Initial Condition h=0.2; % Time step t = 0:h:2; % t goes from 0 to 2 seconds. yexact = 3*exp(-2*t) % Exact solution (in general we won't know this ystar = zeros(size(t)); % Preallocate array (good coding practice) ystar(1) = y0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = -2*ystar(i) % Approx for y gives approx for deriv y1 = ystar(i)+k1*h/2; % Intermediate value (using k1) k2 = -2*y1 % Approx deriv at intermediate value. y2 = ystar(i)+k2*h/2; % Intermediate value (using k2) k3 = -2*y2 % Another approx deriv at intermediate value. y3 = ystar(i)+k3*h; % Endpoint value (using k3) k4 = -2*y3 % Approx deriv at endpoint value. ystar(i+1) = ystar(i) + (k1+2*k2+2*k3+k4)*h/6; % Approx soln end plot(t,yexact,t,ystar); legend('Exact','Approximate');

The MATLAB commands match up easily with the steps of the algorithm. A slight variation of the code was used to show the effect of the size of *h*on the accuracy of the solution (see image below). Note that larger values of *h* result in poorer approximations, but that the solutions are much better than those obtained with the Second Order Runge-Kutta for the same value of *h*. This isn't obvious for small *h*, but is for the curve where *h=0.8*.

For a first order ordinary differential equation defined by

$${{dy(t)} \over {dt}} = f(y(t),t)$$to progress from a point at *t=t₀*, *y*(t₀)*, by one time step, *h*, follow these steps (repetitively).

{k_1} &= f\left( {{y^*}({t_0}),{t_0}} \right)\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{1}}}{\rm{)}} \cr

{k_2} &= f\left( {{y_1}\left( {{t_0} + {h \over 2}} \right),{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y_2}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_2}{h \over 2}\quad &{\rm{another}}\;{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2\;{\rm{(using}}\;{{\rm{k}}_{\rm{2}}}{\rm{)}} \cr

{k_3} &= f\left( {{y_2}\left( {{t_0} + {h \over 2}} \right),{t_0} + {h \over 2}} \right)\quad &{\rm{another}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y_3}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_3}h\quad &{\rm{an}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h\;{\rm{(using}}\;{{\rm{k}}_{\rm{3}}}{\rm{)}} \cr

{{\rm{k}}_4} &= f\left( {{y_3}\left( {{t_0} +h} \right),{t_0} +h} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h \cr \cr

{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {{{k_1} + 2{k_2} + 2{k_3} + {k_4}} \over 6}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$

Notes:

- an initial value of the function must be given to start the algorithm.
- see the MATLAB program on this page for an example.

The global error of the Fourth Order Runge-Kutta algorithm is *O(h ^{4})*.

This is not proven here, but the proof is similar to that for the Second Order Runge-Kutta

The Second Order Runge-Kutta had more than one form (because the technique is derived from an underspecified set of equations). Likewise, the Fourth Order Runge-Kutta has (infinitely many) other forms. A common one is given below, without proof, simply to show that other (possibly very different) forms of the Fourth Order Rung-Kutta can be formulated:

$$\eqalign{{k_1} &= f({y^*}({t_0}),{t_0}) \cr

{k_2} &= f\left( {{y^*}({t_0}) + {k_1}{h \over 3},\;{t_0} + {h \over 3}} \right) \cr

{k_3} &= f\left( {{y^*}({t_0}) - {k_1}{h \over 3} + {k_2}h,\;{t_0} + {2 \over 3}h} \right) \cr

{k_4} &= f\left( {{y^*}({t_0}) + {k_1}h - {k_2}h + {k_3}h,\;{t_0} + h} \right) \cr

{y^*}({t_0} + h) &= {y^*}({t_0}) + {{{k_1} + 3{k_2} + 3{k_3} + {k_4}} \over 8}h \cr} $$

© Copyright 2005 to 2015 Erik Cheever This page may be freely used for educational purposes.

Erik Cheever Department of Engineering Swarthmore College