MATLAB TUTORIAL for the First Course. Part 3: Runge--Kutta 2

Prof. Vladimir A. Dobrushkin

This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Runge--Kutta Methods of order 2
The second order Runge--Kutta method (denoted RK2) simulates the accuracy of the Tylor series method of order 2. Although this method is not as good as the RK4 method, its derivation illustrates all steps and the principles involved. We consider the initial value problem \( y' = f(x,y) , \quad y(x_0 ) = y_0 \) that is assumed to have a unique solution. Using a uniform grid of points \( x_n = x_0 + n\, h , \) with fixed step size h, we write down the Taylor series formula for y(x+h):
\[ y(x+h) = y(x) + h\,y' (x) + \frac{h^2}{2}\, y'' (x) + C\,h^3 + \cdots , \]
where C is a constant expressed through the third derivative of y(x) and the other terms in the series involve powers of h greater than 3.

The derivatives y' and y'' must be extressed in terms of f(x,y) and its partial derivatives according to the given differential equation \( y' = f(x,y) . \) So we have

\begin{align*} y' (x) &= f(x,y) , \\ y'' (x) &= f_x (x,y) + f(x,y)\,f_y (x,y) , \end{align*}
where \( f_x = \partial f/\partial x \quad \mbox{and} \quad f_y = \partial f/\partial y . \) Using these formulas, we derive the Taylor expression for y(x+h):
\[ y(x+h) = y(x) + h\,y' (x) + h\,b_1 f_0 + h\,b_2 f_1 , \]
\[ \begin{split} f_0 &= f(x,y) \\ f_1 &= f(x+c_2 h , y + a_{2,1} h\, f_0 ) . \end{split} \]
Compare the above Runge-Kuta approximation of order 2 with the Taylor polynomial approximation
\begin{align*} y(x+h ) &= y(x) + h\,f(x,y) + \frac{h^2}{2}\, y'' (x) + \cdots , \\ y(x+h ) &= y(x) + \left( b_1 + b_2 \right) h\,f(x,y) + b_2 c_2 h^2 f_x (x,y) + b_2 a_{2,1} h^2 f_y (x,y)\,f(x,y) + b_2 C \,h^3 + \cdots , \end{align*}
we derive the following conclusions:
\begin{align*} h\,f(x,y) &= \left( b_1 + b_2 \right) h\,f(x,y) \qquad \Longrightarrow \qquad &1=b_1 + b_2 , \\ \frac{h^2}{2} \,f_x (x,y) &= h^2 b_2 c_2 f_x (x,y) \qquad \Longrightarrow \qquad & \frac{1}{2} = b_2 c_2 , \\ \frac{h^2}{2} \,f_y (x,y)\, f(x,y) &= h^2 b_2 a_{2,1} f_y (x,y)\, f(x,y) \qquad \Longrightarrow \qquad & \frac{1}{2} = b_2 a_{2,1} . \end{align*}
Hence, if we require that coefficients in the RK2 method satisfy the relations
\[ b_1 + b_2 =1 , \quad \frac{1}{2} = b_2 c_2 , \quad \frac{1}{2} = b_2 a_{2,1} , \]
Then the RK2 method will have the same order of accuracy as the Taylor's method, namely, 2.

Since there are only three equations in four unknowns, the system is underdetermined and we are permitted to choose one of the coefficients as we wish. There are several special choices that have been studied in the literature, and we mention some of them.

If we choose b1 = 1/2, this yields b2 = 1/2, c2 = 1, and a2,1 = 1. So we get Heun's method.

If we choose b1 = 0, this yields b2 = 1, c2 = 1/2, and a2,1 = 1/2. So we get the modified Euler method.

Example. The Butcher table for the explicit midpoint method is given by:

\[ \begin{array}{c|cc} 0&0&0 \\ \frac{1}{2} & \frac{1}{2} & 0 \\ \hline & 0 & 1 \end{array} \]

function  Complete


			Complete


                

 

I. Second Order Optimal Runge-Kutta Approximation


Ralston's method is a second-order method with two stages and a minimum local error bound. Its Mathematica realization is presented below when the step size is denoted by h: \begin{equation} \label{E35.opt} y_{n+1} =y_n + \frac{h}{4}\,f(x_n , y_n ) + \frac{3h}{4} \,f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation} Its Butcher tableau is

\[ \begin{array}{c|cc} 0&0&0 \\ \frac{2}{3} & \frac{2}{3} & 0 \\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} . \]

 

II. Implicit Runge-Kutta of order 2


\begin{equation} \label{E36.implicit} y_{n+1} = y_n + \frac{h}{2}\left( k_1 + k_2 \right) , \end{equation} where \begin{align*} k_1 &= f\left( x_n + \frac{\sqrt{3} -1}{2\sqrt{3}}\,h , \ y_n + \frac{h}{4}\,k_1 + \frac{\sqrt{3} -2}{4\sqrt{3}}\,hk_2\right) , \\ k_2 &= f\left( x_n + \frac{\sqrt{3} +1}{2\sqrt{3}}\,h , \ y_n + \frac{\sqrt{3} +2}{4\sqrt{3}}\,hk_1 + \frac{h}{4}\,k_2\right) , \end{align*}

 

III. Ralston's method


The method was proposed by Anthony Ralston (born in 1930, in New York City). Anthony Ralston received his Ph.D. from the Massachusetts Institute of Technology in 1956 under the direction of M. Eric Reissner. For approximating the solution of the initial value problem
\[ y' = f \left( x, y \right) , \qquad y(x_0 ) = y_0 . \]
Ralston's second order method gives:
\[ y_{i+1} = y_i + \frac{h}{3} \left( k_1 + 2\,k_2 \right) , \qquad x_i = x_0 + i\,h, \]
where
\[ \begin{split} k_1 &= f\left( x_i , y_i \right) , \\ k_2 &= f\left( x_i + \frac{3}{4}\, h , y_i +\frac{3}{4}\, h \, k_1 \right) . \end{split} \]
The Butcher tableau for Ralston's algorithm is
\[ \begin{array}{c|cc} 0&0&0 \\ \frac{3}{4} & \frac{3}{4} & 0 \\ \hline & \frac{1}{3} & \frac{2}{3} \end{array} . \]

The second order Ralson approximation can written as follows:

\begin{equation*} %\label{E35.ralson} y_{n+1} =y_n + \frac{h}{3}\,f(x_n , y_n ) + \frac{2h}{3} \,f\left( x_n + \frac{3h}{4} , y_n + \frac{3h}{4} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation*}

 

IV. Generic second-order method


There is a family of Runge--Kutta methods of order 2, parameterized by an arbitrary positive number α from the interval (0,1] and given by the formula
\[ y_{n+1} = y_n + h \left( 1 - \frac{1}{2\alpha} \right) f(x_n , y_n ) + \frac{1}{2\alpha} \, f \left( x_n + \alpha \,h , y_n + \alpha \,h \,f(x_n , y_n ) \right) , \]
with Butcher tableau
\[ \begin{array}{c|cc} 0&0&0 \\ \alpha & \alpha & 0 \\ \hline & \left( 1 - \frac{1}{2\alpha} \right) & \frac{1}{2\alpha} \end{array} . \]
Correspondingly,
\[ y_{i+1} = y_i + h \left( b_1 k_1 + b_2 k_2 \right) , \qquad b_1 + b_2 =1, \]
where
\[ \begin{split} k_1 &= f\left( x_i , y_i \right) , \\ k_2 &= f\left( x_i + a\,h , y_i +a \, k_1 h \right) . \end{split} \]

This method unites all second-order Runge--Kutta methods; in particular, we have \( a= 1 \) for Heun's algorithm, \( a= 1/2 \) for the midpoint method, and \( a= 2/3 \) for optimal second order Ralston's algorithm.

	function [y, t] = implicit_rk2_ode(f, t0, y0, h, N)
% RK4_ODE Solve a first order diff eq with implicit second order Runge-Kutta method
%   [y, t] = IMPLICIT_RK4_ODE(f, t0, y0, h, N) solves the IVP 
%   {y' = f(t,y), y(t0) = y0} on the interval [t, t + h*N]
%   Example usage:
%   We want to solve the IVP: y' = 1 - t + 4y, y(0) = 1
%   f = @(t,y) 1 - t + 4*y;
%   t0 = 0; y0 = 1; h = 0.001; N = 1000;
%   [y, t] = implicit_rk2_ode(f, t0, y0, h, N);
%   plot(t,y)
%   format long g
%   disp [t' y']

    t = t0:h:(t0+N*h); % set up t-points
    y(N+1)=0; % memory allocation for y-points
    
    % Important note: because Matlab starts with index 1, what we refer to as
    % t_j is called t(j +1) in Matlab; similarly, what we refer to as y_j
    % is called y(j +1) in Matlab
    
    % First y-point is from our initial condition:
    y(1) = y0;
    
    % Now, the main loop of the algorithm:
    for j = 0:N-1
        t_j = t(j +1);
        y_j = y(j +1);
        
        % First, find k1 and k2:
        [k1, k2] = solve_implicit_system(f, h, t_j, y_j);
        
        % Now use the Runge-Kutta equation to find y_j+1:
        y_jplus1 = y_j + (h/2)*(k1+k2);
        
        y(j+2) = y_jplus1;
    end
end	  
Example to run:
f = @(t,y) 1 - t + 4*y;
t0 = 0; y0 = 1; h = 0.001; N = 1000;
[y, t] = implicit_rk2_ode(f, t0, y0, h, N);
plot(t,y)
format long g
disp([t' y'])
To see just one numerical value, type:
 y(t==1.0) 				
Here is another version of implimentation of the second order implicit Runge--Kutta algorithm:
 function rukun
f=inline('-x^2+y') % the right side of ODE
h=0.1;
x0=0; % intial argument
y0=-1; %initial value y(x0)
xn=1; % terminal argument
x=x0:h:xn;
y(numel(x))=0; %memory allocation
y(1)=y0;
for i=1:numel(x)-1 
  k=ss(f, h, x(i), y(i));
  y(i+1)=y(i)+1/2*h*(k(1)+k(2));
end
disp([x',y'])
plot(x,y)
function solution=ss(f,h,x,y)
  s=sqrt(3);
  t1=(s-1)/(2*s);
  t2=(s-2)/(4*s);
  t3=(s+1)/(2*s);
  t4=(s+2)/(4*s);
  F=@(k) [f(x+t1*h, y+(h/4)*k(1)+t2*h*k(2))-k(1);
  f(x+t3*h, y+t4*h*k(1)+h/4*k(2))-k(2)];   
  initial_guess=[1;1];
  [solution]=fsolve(F, initial_guess);