Second Order Runge-Kutta

Contents

Statement of Problem

Consider the situation in which the solution, y(t), to a differential equation

$${{dy(t)} \over {dt}} = y'(t) = f(y(t),t), \quad \quad {\rm{with\;}} y(t_0)=y_0$$

is to be approximated by computer (starting from some known initial condition, y(t₀)=y₀; also, note that the tick mark denotes differentiation). The following text develops an intuitive technique for doing so, and presents some examples. This technique is known as "Second Order Runge-Kutta".

Second Order Runge-Kutta Method (Intuitive)

A First Order Linear Differential Equation with No Input

The first order Runge-Kutta method used the derivative at time t₀ (t₀=0 in the graph below) to estimate the value of the function at one time step in the future. If you are not familiar with it, you should read the section entitled: A First Order Linear Differential Equation with No Input. We repeat the central concept of generating a step forward in time in the following text.

$${{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. The exact solution in this case is y(t)=3e-2t, t≥0, though in general we won't know this and will need numerical integration methods to generate an approximation.

In the graph below, the slope at t=0 is called k1, and the estimate is called y*(h); in this example h=0.2.

This obviously leads to some error in the estimate, and we would like to reduce this error. One way we could do this, conceptually, is to use the derivative at the halfway point between t=0 and t=h=0.2. The slope at this point (t=½h=0.1) is shown below (and is labeled k2). Note the line (orange) is tangent to the curve (blue) at t=½h.

Now if we use this intermediate slope, k2, as we step ahead in time then we get better estimate, y*(h), than we did before. On the diagram below the exact value of the solution is y(0.1)=2.0110 and the approximation is y*(0.1)=2.0175 for an error of about 0.3% (compared with about 10% error for the first order Runge-Kutta).

This seems like a very nice solution, and obviously generates a significantly more accurate approximation than the first order technique that uses a line with slope, k1, calculated at t=0. The problem is we don't know the exact value of y(½h) so we can't find the exact value of k2the slope at t=½h (Recall that the calculation of the derivative requires knowledge of the value of the function, y'(t) = - 2y(t)).

What we do instead is use the First Order Runge-Kutta to generate an approximate value for y(t) at t=½h=0.1, call it y1(½h). We then use this estimate to generate k2 (which will be an approximation to the slope at the midpoint), and then use k2 to find y*(h). To step from the starting point at t=0 to an estimate at t=h, follow the procedure below.

$$\eqalign{
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 \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(h) &= y(0) + y'(0)h + y''(0){{{h^2}} \over 2} + \cdots \quad &{\rm{Taylor}}\;{\rm{Series}}\;{\rm{around}}\;t = 0 \cr
y(t) &\approx y(0) + y'(0)h\quad &{\rm{Truncate}}\;{\rm{Taylor}}\;{\rm{Series}} \cr
{y^*}(h) &= y(0) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y(h) \cr} $$

In general, to go from the estimate t=t₀ to an estimate at t=t₀+h

$$\eqalign{
y'({t_0}) &= - 2y({t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{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 \cr
{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr
y({t_0} + h) &= y({t_0}) + y'({t_0})h + y''(0){{{h^2}} \over 2} + \cdots \quad &{\rm{Taylor}}\;{\rm{Series}}\;{\rm{around}}\;t = {t_0} \cr
y({t_0} + h) &\approx y({t_0}) + y'({t_0})h\quad &{\rm{Truncated}}\;{\rm{Taylor}}\;{\rm{Series}} \cr
{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y(t_0+h) \cr} $$
Example 1: Approximation of First Order Differential Equation with No Input Using MATLAB

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

%% Example 1
% Solve y'(t)=-2y(t) with y0=3, midpoint method
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
  k2 = -2*y1                    % Approx deriv at intermediate value.
  ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
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 hon 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 First Order Runge-Kutta for the same value of h.

Key Concept: Second Order Runge-Kutta Algorithm (midpoint)

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).

$$\eqalign{
{k_1} &= f({y^*}({t_0}),\;{t_0})\quad &{\rm{estimate}}\;{\rm{of}}\;{\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 \over 2} \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 \over 2} \cr
{y^*}\left( {{t_0} + h} \right)& = {y^*}({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y\left( {{t_0} + h} \right) \cr} $$

Notes:

A First Order Linear Differential Equation with Input

Adding an input function to the differential equation presents no real difficulty. You just need to ensure that you evaluate the function at t=t₀ to find k1, and at t=t₀+frac12;h to find k2. Consider an input of cos(4t).

$$\eqalign{ y'(t) + 2y(t) &= \cos (4t) \quad {\rm{or}} \quad y'(t) &= - 2y(t) + \cos (4t) \cr} $$

So

$$\eqalign{
y'({t_0}) &= - 2y({t_0}) + \cos (4{t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{k_1} &= - 2{y^*}({t_0}) + \cos (4{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 \cr
{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right) + \cos \left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr
{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$

Note: the exact solution in this case is $y(t) = 2.9{e^{ - 2{\kern 1pt} t}} + 0.1\cos (4t) + 0.2\sin (4t)$.

Example 2: Approximation of First Order Differential Equation with Input Using MATLAB

We can use MATLAB to perform the calculation described above. To perform this new approximation all that is necessary is to change the calculation of k1 and k2 using the appropriate value for the time variable (the value of the exact solution is also changed, for plotting).

%% Example 2
% Solve y'(t)=-2y(t)-cos(4t) with y0=3
y0 = 3;                  % Initial Condition
h=0.2;                   % Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
yexact = 0.1*cos(4*t)+0.2*sin(4*t)+2.9*exp(-2*t);  % Exact Solution     
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)+cos(4*t(i));  % Approx for y gives approx for deriv
  y1 = ystar(i)+k1*h/2;      % Intermediate value
  k2 = -2*y1+cos(4*(t(i)+h/2));        % Approx deriv at intermediate value.
  ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A modified version of the program (that uses several values of h) generated the plot below. As before, the solution is better with smaller values of h, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of h.

A First Order Non-Linear Differential Equation

Using a non-linear differential equation presents no real difficulty. Again, you just need to ensure that you evaluate the function at t=t₀ to find k1, and at t=t₀+frac12;h to find k2. Consider

$${{dy(t)} \over {dt}} - y(t)(1 - 2t) = 0,\quad \quad \quad \quad y(0) = 1$$

with solution (described here)

$$y(t) = {e^{t - {t^2}}}$$

The process follows those that came before

$$\eqalign{
y'({t_0}) &= y({t_0})(1 - 2{t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{k_1} &= y^*({t_0})(1 - 2{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 \cr
{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right)\left( {1 - 2\left( {{t_0} + {h \over 2}} \right)} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr
{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$
Example 3: Approximation of First Order Nonlinear Differential Equation with Input Using MATLAB

As before, to perform this new approximation all that is necessary is to change the calculation of k1 and k2 using the appropriate value for the time variable (the value of the exact solution is also changed, for plotting).

%% Example 3
% Solve y'(t)=y(t)(1-2t) with y0=1
y0 = 1;                  % Initial Condition
h=0.2;                   % Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
% Exact solution, hard to find in this case (in general we won't have it)
yexact = exp(t-t.^2);
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 = ystar(i)*(1-2*t(i));     % Approx for y gives approx for deriv
  y1 = ystar(i)+k1*h/2;         % Intermediate value
  k2 = y1*(1-2*(t(i)+h/2));     % Approx deriv at intermediate value.
  ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A modified version of the program (that uses several values of h) generated the plot below. As before, the solution is better with smaller values of h, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of h.

A Higher Order Linear Differential Equation

If we start with a higher order differential equation

$$\displaylines{
{{{d^3}y(t)} \over {dt}} + 4{{{d^2}y(t)} \over {dt}} + 6{{dy(t)} \over {dt}} + 4y(t) = \gamma (t)\quad \quad \quad \gamma (t) = {\rm{unit step function}} \cr
{\left. {{{{d^2}y(t)} \over {dt}}} \right|_{t = {0^ + }}} = 0,\quad \quad {\left. {{{dy(t)} \over {dt}}} \right|_{t = {0^ + }}} = - 1,\quad \quad y({0^ + }) = 0\quad \cr} $$

We can rewrite as a matrix differential equation (see previous example if this is unclear)

$$\eqalign{ {\bf{q}}'(t) &= \left[ {\matrix{ {{q_1}'(t)} \cr {{q_2}'(t)} \cr {{q_3}'(t)} \cr} } \right] = \left[ {\matrix{ 0 & 1 & 0 \cr 0 & 0 & 1 \cr
{ - 4} & { - 6} & { - 4} \cr} } \right]\left[ {\matrix{ {{q_1}(t)} \cr {{q_2}(t)} \cr {{q_3}(t)} \cr} } \right] + \left[ {\matrix{ 0 \cr 0 \cr 1 \cr} } \right]\gamma (t) \cr
{\bf{q}}'(t) &= {\bf{Aq}}(t) + {\bf{B\gamma }}(t) \cr {\bf{A}} &= \left[ {\matrix{ 0 & 1 & 0 \cr 0 & 0 & 1 \cr
{ - 4} & { - 6} & { - 4} \cr} } \right],\quad \quad \quad {\bf{B}} = \left[ {\matrix{ 0 \cr 0 \cr 1 \cr } } \right] \cr} $$

Proceed as before (using matrices in place of scalars where appropriate

$$\eqalign{
{\bf{q}}'({t_0}) &= {\bf{Aq}}({t_o}) + {\bf{B}}\gamma ({t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{{\bf{k}}_1} &= {\bf{A}}{{\bf{q}}^*}({t_o}) + {\bf{B}}\gamma ({t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{{\bf{q}}_1}\left( {{t_0} + {h \over 2}} \right) &= {{\bf{q}}^*}({t_o}) + {{\bf{k}}_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2 \cr
{{\bf{k}}_2} &= {\bf{A}}{{\bf{q}}_1}\left( {{t_0} + {h \over 2}} \right) + {\bf{B}}\gamma \left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr
{{\bf{q}}^*}\left( {{t_0} + h} \right) &= {{\bf{q}}^*}({t_o}) + {{\bf{k}}_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;{\bf{q}}({t_0} + h) \cr} $$

Note: the exact solution to this problem is $y(t) = {1 \over 4}+{{\rm{e}}^{ - t}}{\mkern 1mu} \left( {\cos (t) - {5 \over 2}\sin (t)} \right) - {5 \over 4}{{\rm{e}}^{ - 2{\kern 1pt} t}}$, for t≥0.

Example 4: Approximation of Third Order Differential Equation Using MATLAB

To perform this new approximation all that is necessary is to change the approriate variables from scalars to vectors or matrices, and to define the A and B matrices. Note that in this case we multiply the B matrix by 1 since the input is a unit step (γ(t)=1 for t≥0). If the input were a sine wave, the sine wave would multiply B; if there is no input, it is not necessary to include the B matrix at all (it is multiplied by 0).

%% Example 4
% Solve y'''(t)+4y''(t)+6y'(t)+4y(t)=gamma(t),  y''(0)=0, y'(0)=-1, y(0)=0
q0 = [0; -1; 0];         % Initial Condition (vector)
h=0.2;                   % Time step
t = 0:h:5;               % t goes from 0 to 2 seconds.
A = [0 1 0; 0 0 1; -4 -6 -4];  % A Matrix
B = [0; 0; 1];                 % B Matrix

yexact = 1/4 + exp(-t).*(cos(t) - 5*sin(t)/2) - 5*exp(-2*t)/4;  % Exact soln
qstar = zeros(3,length(t));  % Preallocate array (good coding practice)

qstar(:,1) = q0;             % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
  k1 = A*qstar(:,i)+B*1;            % Approx for y gives approx for deriv
  q1 = qstar(:,i)+k1*h/2;           % Intermediate value
  k2 = A*q1+B*1;                    % Approx deriv at intermediate value.
  qstar(:,i+1) = qstar(:,i) + k2*h; % Approx solution at next value of q
end
plot(t,yexact,t,qstar(1,:));        % ystar = first row of qstar
legend('Exact','Approximate');

A modified version of the program (that uses several values of h) generated the plot below. As expected, the solution is better with smaller values of h, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of h.

Second Order Runge-Kutta Method (The Math)

The Second Order Runge-Kutta algorithm described above was developed in a purely ad-hoc way. It seemed reasonable that using an estimate for the derivative at the midpoint of the interval between t₀ and t₀+h (i.e., at t₀+½h) would result in a better approximation for the function at t₀+h, than would using the derivative at t₀ (i.e., Euler's Method &emdash; the First Order Runge-Kutta). And, in fact, we showed that the resulting approximation was better. But, is there a way to derive the Second Order Runge-Kutta from first principles? If so, we might be able to develop even better algorithms.

Aside: Some Useful Math Review

In the following derivation we will use two math facts that are reviewed here. You should be familiar with this from a course in multivariable calculus.


1) If there is a function of two variable, g(x,y), it's Taylor Series about x₀, y₀ is given by:

$$\eqalign{ g({x_0} + \Delta x,{y_0} + \Delta y) &= g({x_0},{y_0}) + {\left. {{{\partial g(x,y)} \over {\partial x}}} \right|_{{x_o},{y_0}}}\Delta x + {\left. {{{\partial g(x,y)} \over {\partial y}}} \right|_{{x_o},{y_0}}}\Delta y + \cdots \cr
&= g + {g_x}\Delta x + {g_y}\Delta y + \cdots \cr} $$

where Δx is the increment in the first dimenstion, and Δy is the increment in the second dimension. In the last line we use a shorthand notation that removes the explicit functional notation, and also represents the partial derivative of g with respect to x as gx (and likewise for gy).


2) If z is a function of two variables z=g(x,y), where x=r(t) and y=s(t), then by the chain rule for partial derivatives

$${{dz} \over {dt}} = {{\partial g} \over {\partial r}}{{dr} \over {dt}} + {{\partial g} \over {\partial s}}{{ds} \over {dt}}$$

In particular if

$$y'(t) = {{dy(t)} \over {dt}} = f(y(t),t) = f$$

then

$$y''(t) = {{df(y(t),t)} \over {dt}} = {{\partial f} \over {\partial t}}{{dt} \over {dt}} + {{\partial f} \over {\partial y}}{{dy} \over {dt}} = {{\partial f} \over {\partial t}} + {{\partial f} \over {\partial y}}{{dy} \over {dt}} = {f_t} + {f_y}f$$

To start, recall that we are expressing our differential equation as

$${{dy(t)} \over {dt}} = y'\left( t \right) = f(y(t),t)$$

Now define two approximations to the derivative

$$\eqalign{
& {k_1} = f(y^*({t_0}),{t_0}) \cr
& {k_2} = f(y^*({t_0}) + \beta h{k_1},{t_0} + \alpha h) \cr} $$

In all cases α and β will represent fractional values between 0 and 1. These equation state that k1 is the approximation to the derivative based on the estimated value of y(t) at t=t₀ (i.e., y*(t₀)) and the time at t₀. The value of k2 is based upon the estimated value, y*(t₀), plus some fraction of the step size, βh, times the slope k1, and the time at t₀+αh (i.e., some time between t₀ and t₀+h). In the method described previously α=β=½, but other choices are possible.

To update our solution with the next estimate of y(t) at t₀+h we use the equation

$$y^*({t_0} + h) = y^*({t_0}) + h(a{k_1} + b{k_2})\quad \quad \quad {\rm{update\;equation}}$$

This equation states that we get the value of y*(t₀+h) from the value of y*(t₀) plus the time step, h, multiplied by a slope that is a weighted sum of k1 and k2. In the method described previously a=0 and b=1, so we used only the second estimate for the slope. (Note that Euler's Method (First Order Runge-Kutta) is a special case of this method with a=1, b=0, and α and β don't matter because k2 is not used in the update equation.)

Our goal now is to determine, from first principles, how to find the values a, b, α and β that result in low error. Starting with the update equation (above) we can substitue the previously given expressions for k1 and k2 which yields

$${y^*}({t_0} + h) = {y^*}({t_0}) + h\left( {af\left( {{y^*}({t_0}),{t_0}} \right) + bf\left( {{y^*}({t_0}) + \beta h{k_1},\,{t_0} + \alpha h} \right)} \right)$$

We can use a two-dimensional Taylor Series (where the increment in the first dimension is βhk1, and the increment in the second dimension is αh) to expand the rightmost term.

$$\eqalign{
f\left( {{y^*}({t_0}) + \beta h{k_1},\,{t_0} + \alpha h} \right) &= f + {f_y}\beta h{k_1} + {f_t}\alpha h + \cdots \cr
&= f + {f_y}f\beta h + {f_t}\alpha h + \cdots \cr} $$

In the last line we used the fact the k1=f . The ellipsis denotes terms that are second order or higher. Substituting this into the previous expression for y*(t₀+h) and rearranging we get

$$\eqalign{
{y^*}({t_0} + h) &= {y^*}({t_0}) + h\left( {af + b\left( {f + {f_y}f\beta h + {f_t}\alpha h + \cdots } \right)} \right) \cr
&= {y^*}({t_0}) + haf + hbf + {h^2}b\beta {f_y}f + {h^2}{f_t}b\alpha + \cdots \cr
&= {y^*}({t_0}) + h(a + b)f + {h^2}{f_t}b\alpha + {h^2}b\beta {f_y}f + \cdots \cr} $$

The ellipsis was multiplied by h between the first or second line and now represents terms that are third order or higher.

To finish we compare this approximation with the expressionfor a Taylor Expansion of the exact solution (going from the first line to the second we used the chain rule for partical derivatives). In this expression the ellipsis represents terms that are third order or higher.

$$\eqalign{
y({t_0} + h) &= y({t_0}) + h{\left. {y'(t)} \right|_{{t_0}}} + {{{h^2}} \over 2}{\left. {y''(t)} \right|_{{t_0}}} + \cdots \cr
&= y({t_0}) + hf + {{{h^2}} \over 2}\left( {{f_t} + {f_y}f} \right) + \cdots \cr
&= y({t_0}) + hf + {{{h^2}} \over 2}{f_t} + {{{h^2}} \over 2}{f_y}f + \cdots \cr} $$

Comparing this expression with our final expression for the approximation,

$${y^*}({t_0} + h) = {y^*}({t_0}) + hf(a + b) + {h^2}{f_t}b\alpha + {h^2}b\beta {f_y}f + \cdots $$

we see that they agree up to the error terms (third order and higher) represented by the ellipsis if we define the constants, a, b, α and β such that

$$\displaylines{
a + b = 1 \cr
b\alpha = {1 \over 2} \cr
b\beta = {1 \over 2} \cr} $$

This system is underspecified, there are four unknowns, and only three equations, so more than one solution is possible. Clearly the choice we made a=0 and b=1 and α=β=½ is one of these solutions. Because all of the terms of the approximation are equal to the terms in the exact solution, up to the error terms, the local error of this method is therefore O(h3) (O(h2) globally, hence the term "second order" Runge-Kutta).

Key Concept: Error of Second Order Runge Kutta

The global error of the Second Order Runge-Kutta algorithm is O(h2).

Another Form of the Second Order Runge-Kutta Method

Another common choice for the coefficients of the algorithm are a=b=½ and α=β=1. Before giving an example, let's figure out, intuitively what this is doing. We start with our equations for k1, k2, and y*(t₀+h).

$$\eqalign{
{k_1} &= f({y^*}({t_0}),{t_0}) \cr
{k_2} &= f({y^*}({t_0}) + \beta h{k_1},{t_0} + \alpha h) \cr
{y^*}({t_0} + h) &= {y^*}({t_0}) + h(a{k_1} + b{k_2}) \cr} $$

Now put in our chosen constants

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

In terms of our algorithm description

$$\eqalign{
{k_1} &= f({y^*}({t_0}),\;{t_0})\quad& {\rm{estimate}}\;{\rm{of}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{y_1}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_1}h\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h \cr
{k_2} &= f\left( {{y_1}\left( {{t_0} + h} \right),\;{t_0} + h} \right)\quad& {\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h \cr
{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + h\left( {{{{k_1} + {k_2}} \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{y}}\left( {{{\rm{t}}_{\rm{0}}}{\rm{ + h}}} \right)\;{\rm{using}}\;{\rm{average}}\;{\rm{of}}\;{\rm{slopes}} \cr} $$

We will again use the 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 with exact solution y(t)=3e-2t, t≥0. Conceptually what we are trying to do is to find the derivative at the the beginning of our time step (t=0, for the first step) and at the end of our time step, t=h=0.2. We call the slope at the beginning k1; numerically k1=-6. The slope at this point is shown below (and is labeled k1). Note the line (orange) is tangent to the curve (blue) at t=0.

The slope at the end point is shown below (and is labeled k2); numerically k2=-4.0219. Note the line (orange) is tangent to the curve (blue) at t=h; in this example h=0.2.

We now take the average of these two slopes (-5.0110) and use it to move forward by one time step.

$${y^*}\left( h \right) = {y^*}(0) + h\left( {{{{k_1} + {k_2}} \over 2}} \right) = 3 + 0.2\left( { - 5.0110} \right) = 1.9978$$

This is very close to the exact answer, y(0.2)=2.0110, with an error of about 0.6%. This is close to that of the midpoint method (it is slightly worse, but that is becaus of the specific problem chosen, that won't generally be the case), and much better than that of the first order algorithm.

Key Concept: The Second Order Runge-Kutta Algorithm (endpoint)

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).

$$\eqalign{
{k_1} &= f({y^*}({t_0}),\;{t_0})\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr
{y_1}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_1}h\quad &{\rm{approximate}}\;{\rm{enpoint}}\;{\rm{value}}\;{\rm{at}}\;t = {t_0} + h \cr
{k_2} &= f\left( {{y_1}\left( {{t_0} + h} \right),\;{t_0} + h} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;{\rm{enpoint,}}\;t = {t_0} + h \cr
{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + h{{\left( {{k_1} + {k_2}} \right)} \over 2}\quad &{\rm{estimate}}\;{\rm{of}}\;y\left( {{t_0} + h} \right)\;{\rm{using}}\,{\rm{average}}\;{\rm{slope}} \cr} $$

 

Example 5: Repeat Example 1 with endpoints method

The following MATLAB code repeats Example 1 (a linear differential equation with no input). Example 1 used the "midpoint" method, this example uses the "endpoint" method. The MATLAB commands match up easily with the steps of the algorithm (only the lines that calculate y1 and k2 have changed from the midpoint method).

%% Example 5
% Solve y'(t)=-2y(t) with y0=3, endpoint method
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;           % enpoint approximation 
  k2 = -2*y1                    % Approx deriv at endpoint.
  ystar(i+1) = ystar(i) + h*(k1+k2)/2; % Approx solution at next value of y
end
plot(t,yexact,t,ystar);
legend('Exact','Approximate');

A slight variation of the code was used to show the effect of the size of hon 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 First Order Runge-Kutta for the same value of h (and appears similar to results obtained with the mid-point method).

The other examples can be coded for the endpoint method in a similar way.

Moving on

With a little more work we can develop an algorithm that is accurate to even higher order. The next page describes just such a method.


References

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

Erik Cheever       Department of Engineering         Swarthmore College