http://lpsa.swarthmore.edu/NumInt/NumIntFourth.html

Fourth Order Runge-Kutta

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

d y ( t ) d t = y' (t) = f (y (t), t), w i t h 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) k 2 = f (y * (t 0) + k 1 h 2, t 0 + h 2) k 3 = f (y * (t 0) + k 2 h 2, t 0 + h 2) k 4 = f (y * (t 0) + k 3 h, t 0 + h)

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)

y * (t 0 + h) = y * (t 0) + k 1 + 2 k 2 + 2 k 3 + k 4 6 h = y * (t 0) + (1 6 k 1 + 1 3 k 2 + 1 3 k 3 + 1 6 k 4) h = y * (t 0) + m h w h e r e m i s a w e i g h t e d a v e r a g e s l o p e a p p r o x i m a t i o n

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

d y ( t ) d t + 2 y (t) = 0 o r d y ( t ) d t = - 2 y (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) k 1 y 1 (h 2) k 2 y 2 (h 2) k 3 y 3 (h) k 4 y * (h) = - 2 y (0) = - 2 y (0) = y * (0) + k 1 h 2 = - 2 y 1 (h 2) = y * (0) + k 2 h 2 = - 2 y 2 (h 2) = y * (0) + k 3 h = - 2 y 3 (h) = y (0) + k 1 + 2 k 2 + 2 k 3 + k 4 6 h e x p r e s s i o n f o r d e r i v a t i v e a t t = 0 d e r i v a t i v e a t t = 0 i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = h / 2 (u s i n g k 1) e s t i m a t e o f s l o p e a t t = h / 2 a n o t h e r i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = h / 2 (u s i n g k 2) a n o t h e r e s t i m a t e o f s l o p e a t t = h / 2 a n e s t i m a t e o f f u n c t i o n a t t = h (u s i n g k 3) e s t i m a t e o f s l o p e a t t = h e s t i m a t e o f y (h)

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

k 1 y 1 (t 0 + h 2) k 2 y 2 (t 0 + h 2) k 3 y 3 (t 0 + h) k 4 y * (t 0 + h) = - 2 y * (t 0) = y * (t 0) + k 1 h 2 = - 2 y 1 (h 2) = y * (t 0) + k 2 h 2 = - 2 y 2 (h 2) = y * (t 0) + k 3 h = - 2 y 3 (t 0 + h) = y * (t 0) + k 1 + 2 k 2 + 2 k 3 + k 4 6 h a p p r o x i m a t e d e r i v a t i v e a t t = t 0 i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = t 0 + h / 2 (u s i n g k 1) e s t i m a t e o f s l o p e a t t = t 0 + h / 2 a n o t h e r i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = t 0 + h / 2 (u s i n g k 2) a n o t h e r e s t i m a t e o f s l o p e a t t = t 0 + h / 2 a n e s t i m a t e o f f u n c t i o n a t t = t 0 + h (u s i n g k 3) e s t i m a t e o f s l o p e a t t = t 0 + h e s t i m a t e o f y (t 0 + 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

d y ( t ) d t = y' (t) = - 2 y (t), w i t h 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₁.

k 1 = - 2 y (t 0) = - 2 y (0) = - 6

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.

k 2 = - 2 y 1 = - 4.8

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

k 3 = - 2 y 2 = - 5.04

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.

k 4 = - 2 y 3 = - 3.9840

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

m = 1 6 k 1 + 1 3 k 2 + 1 3 k 3 + 1 6 k 4

and use it to generate our final estimate

y * (h) = y (0) + (1 6 k 1 + 1 3 k 2 + 1 3 k 3 + 1 6 k 4) 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.

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:

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

Key Concept: Fourth Order Runge-Kutta Algorithm

For a first order ordinary differential equation defined by

d y ( t ) d t = 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 y 1 (t 0 + h 2) k 2 y 2 (t 0 + h 2) k 3 y 3 (t 0 + h) k 4 y * (t 0 + h) = f (y * (t 0), t 0) = y * (t 0) + k 1 h 2 = f (y 1 (t 0 + h 2), t 0 + h 2) = y * (t 0) + k 2 h 2 = f (y 2 (t 0 + h 2), t 0 + h 2) = y * (t 0) + k 3 h = f (y 3 (t 0 + h), t 0 + h) = y * (t 0) + k 1 + 2 k 2 + 2 k 3 + k 4 6 h a p p r o x i m a t e d e r i v a t i v e a t t = t 0 i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = t 0 + h / 2 (u s i n g k 1) e s t i m a t e o f s l o p e a t t = t 0 + h / 2 a n o t h e r i n t e r m e d i a t e e s t i m a t e o f f u n c t i o n a t t = t 0 + h / 2 (u s i n g k 2) a n o t h e r e s t i m a t e o f s l o p e a t t = t 0 + h / 2 a n e s t i m a t e o f f u n c t i o n a t t = t 0 + h (u s i n g k 3) e s t i m a t e o f s l o p e a t t = t 0 + h e s t i m a t e o f y (t 0 + h)

Notes:

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

Key Concept: Accuracy of the Fourth Order Runge-Kutta Algorithm

The global error of the Fourth Order Runge-Kutta algorithm is O(h⁴).
This is not proven here, but the proof is similar to that for the Second Order Runge-Kutta

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:

k 1 k 2 k 3 k 4 y * (t 0 + h) = f (y * (t 0), t 0) = f (y * (t 0) + k 1 h 3, t 0 + h 3) = f (y * (t 0) - k 1 h 3 + k 2 h, t 0 + 2 3 h) = f (y * (t 0) + k 1 h - k 2 h + k 3 h, t 0 + h) = y * (t 0) + k 1 + 3 k 2 + 3 k 3 + k 4 8 h

Diseño Electrónico

sábado, 4 de noviembre de 2017

RK4 Fourth Order Runge-Kutta

Fourth Order Runge-Kutta

Contents

Introduction

The Fourth Order-Runge Kutta Method.