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

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 3

We consider the initial value problem \( y' = f(x,y) , \quad y(x_0 ) = y_0 \) that is assumed to have a unique solution. As usual, we use a uniform grid of points \( x_n = x_0 + n\, h , \quad n= 0,1,2,\ldots , \) with fixed step size h to approximate the actual solution \( y= \phi (x) . \)

We start with the following third order algorithm credited to Kutta:

\begin{equation} \label{E35.9} y_{n+1} = y_n + \frac{h}{6}\, ( k_1 + 4 k_2 + k_3 ), \end{equation} where \[ \begin{array}{l} k_1 = f( x_n , y_n ) , \\ k_2 = f( x_n + h/2 , y_n + k_1 h/2 ) , \\ k_3 = f( x_n +h , y_n + 2h k_2 - h k_1 ) . \end{array} \]

 

The Nystrom Algorithm


The Nystrom algorithm:

\begin{equation} \label{E35.nystrom} y_{n+1} =y_n + \frac{h}{8}\,(2k_1 + 3k_2 + 3k_3 ), \qquad n=0,1,2,\ldots , \end{equation} where \[ k_1 = f(x_n , y_n ), \ k_2 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_1 \right) , \ k_3 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_2 \right) ; \]

 

The Nearly Optimal Algorithm


The nearly optimal algorithm:

\begin{equation} \label{E35.optimal} y_{n+1} =y_n + \frac{h}{9} \,(2k_1 + 3k_2 + 4k_3 ), \qquad n=0,1,2,\ldots , \end{equation} in which the stages are computed according to \[ k_1 = f(x_n , y_n ), \quad k_2 = f \left( x_n + \frac{h}{2} , y_n + \frac{h}{2}\,k_1 \right) , \quad k_3 = f \left( x_n + \frac{3h}{4} , y_n + \frac{3h}{4}\,k_2 \right) ; \]

 

The Heun Algorithm


Heun's algorithm (1900):

\begin{equation} \label{E35.heun} y_{n+1} =y_n + \frac{h}{4} \,(k_1 + 3k_3 ), \qquad n=0,1,2,\ldots , \end{equation}
where
\[ k_1 = f(x_n , y_n ), \ k_2 = f\left( x_n + \frac{h}{3} , y_n + \frac{h}{3}\,k_1 \right) , \ k_3 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_2 \right) . \]

function  Complete


			Complete

 

Bogacki and Shampine Approximation


The Bogacki--Shampine algorithm is a method for the numerical solution of ordinary differential equations, that was proposed in 1989 by Lawrence F. Shampine and his student Przemyslaw Bogacki. Przemyslaw is now a Professor at Old Dominion University. The method is utilized in matlab under the name of ode23. The " 23 " in the function name indicates that two simultaneous single-step formulas, one of second order and one of third order, are involved. The Bogacki–Shampine algorithm is a Runge–Kutta method of order three with four stages, so that it uses approximately three function evaluations per step. It has an embedded second-order method which can be used to implement adaptive step size similar to RKF45. Because of its simplicity, ode23 chooses frequently the natural step size that is well suited for display purposes.

For the initial value problem \( y' = f(x,y) , \quad y(x_0 ) = y_0 , \) the Bogacki–Shampine algorithm determines approximate values yk of the actual solution at mesh points \( x_n = x_0 + n\, h , \quad n=0,1,2,\ldots ; \) as follows:
\begin{align*} k_1 &= f\left( x_n , y_n \right) , \\ k_2 &= f\left( x_n+ \frac{h}{2} , y_n + \frac{h}{2}\, k_1 \right) , \\ k_3 &= f\left( x_n+ \frac{3\,h}{4} , y_n + \frac{3\,h}{4}\, k_2 \right) , \\ y_{n+1} &= y_n + \frac{2}{9}\, h\, k_1 + \frac{1}{3}\,h\,k_2 + \frac{4}{9}\, h\, k_3 , \\ k_4 &= f\left( x_n +h , y_{n+1} \right) , \\ z_{n+1} &= y_n + \frac{7}{24}\,h \,k_1 + \frac{1}{4}\,h\, k_2 + \frac{1}{3}\, h\,k_3 + \frac{1}{8}\, h\, k_4 . \end{align*}
Here zn+1 is a second-order approximation to the exact solution. The third order method for calculating yn+1 is due to Anthony Ralston (1965). He was 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.

The Butcher tableau for the Bogacki–Shampine method is:
\[ \begin{array}{c|cccc} 0&0&0 \\ \frac{1}{2} & \frac{1}{2} & 0 \\ \hline \frac{3}{4} & 0 & \frac{1}{4} && \\ 1 & \frac{2}{9} & \frac{1}{3} & \frac{4}{9} \\ \hline & \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & 0 \\ & \frac{7}{24} & \frac{1}{4} & \frac{1}{3} & \frac{1}{8} \end{array} . \]

References:

Bogacki, Przemyslaw; Shampine, Lawrence F. (1989), "A 3(2) pair of Runge–Kutta formulas", Applied Mathematics Letters, 2 (4): 321--325, ISSN 0893-9659, doi:10.1016/0893-9659(89)90079-7.
Ralston, Anthony (1965), A First Course in Numerical Analysis, New York: McGraw-Hill.
Shampine, Lawrence F.; Reichelt, Mark W. (1997), "The Matlab ODE Suite", SIAM Journal on Scientific Computing, 18 (1): 1–22, ISSN 1064--8275, doi:10.1137/S1064827594276424.