Fourth Order Runge-Kutta

• Intro
• First Order
• Second
• Fourth
• Printable

Contents

Introduction

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)

The Fourth Order-Runge Kutta Method.

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

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

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 k₁ in the first and second order methods).
• If we use the slope k₁ to step halfway through the time step, then k₂ is an estimate of the slope at the midpoint. This is the same as the slope, k₂, from the second order midpoint method. This slope proved to be more accurate than k₁ for making new approximations for y(t).
• If we use the slope k₂ to step halfway through the time step, then k₃ 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 t₀+h), and k₄ is an estimate of the slope at the endpoint.

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

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₂ and k₃) is higher than those at the endpoints (k₁ and k₄), since we expect these to be a better estimate of the slope when going from y*(t₀) to y*(t₀+h).

As an example consider the differential equation

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

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

Visualizing the Fourth Order Runge-Kutta Method

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.

The first slope, k₁ (and finding y₁)

Given y(0), we find the slope, y'(0) using the equation stated above. We call this slope k₁.

We use this slope to estimate the value of the function midway through the intervale, i.e., y(½h). We call this estimate y₁, y₁=y(0)+k₁½h=2.4

The second slope, k₂ (and finding y₂)

Given y₁, 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₂, with y₂=y(0)+k₂½h=2.52

The third slope, k₃ (and finding y₃)

Given y₂, we can use it to find another estimate the slope at t=½h that we will call k₃.

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₃, with y₃=y(0)+k₃h=1.992

The fourth slope, k₄

We use y₃ to find an estimate of the slope at the end of the interval, t=h.

The final slope (a weighted average of previous slopes) — (and finding y*(t₀+h))

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

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.

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

An Alternate form of the Fourth 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:

References