Second Order Runge-Kutta

• Intro
• First Order
• Second
• Fourth
• Printable

Contents

Statement of Problem

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

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.

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 k₁, 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 k₂). Note the line (orange) is tangent to the curve (blue) at t=½h.

Now if we use this intermediate slope, k₂, 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, k₁, 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 k₂the 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 y₁(½h). We then use this estimate to generate k₂ (which will be an approximation to the slope at the midpoint), and then use k₂ to find y*(h). To step from the starting point at t=0 to an estimate at t=h, follow the procedure below.

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

%% 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');

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 k₁, and at t=t₀+frac12;h to find k₂. Consider an input of cos(4t).

So

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
% 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 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 k₁, and at t=t₀+frac12;h to find k₂. Consider

with solution (described here)

The process follows those that came before

%% 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 Higher Order Linear Differential Equation

If we start with a higher order differential equation

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

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

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
% 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');

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.

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

Now define two approximations to the derivative

In all cases α and β will represent fractional values between 0 and 1. These equation state that k₁ 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 k₂ is based upon the estimated value, y*(t₀), plus some fraction of the step size, βh, times the slope k₁, 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

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 k₁ and k₂. 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 k₂ 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 k₁ and k₂ which yields

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

In the last line we used the fact the k₁=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

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.

Comparing this expression with our final expression for the approximation,

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

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(h³) (O(h²) globally, hence the term "second order" Runge-Kutta).

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 k₁, k₂, and y*(t₀+h).

Now put in our chosen constants

In terms of our algorithm description

We will again use the the differential equation

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 k₁; numerically k₁=-6. The slope at this point is shown below (and is labeled k₁). 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 k₂); numerically k₂=-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.

%% 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');

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