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
Each Runge--Kutta method is derived from an appropriate Taylor method in such a way that the final global error is of order O(hm), so it is of order m. These methods can be constructed for any order m. A trade-off to perform several slope evaluations at each step and obtain a good accuracy makes the Runge--Kutta of order 4 the most popular.
As usual, we consider an initial value problem \( y' = f(x,y) , \quad y(x_0 )= y_0 \) and seek approximate values yk to the actual solution \( y= \phi (x) \) at mesh points \( x_n = x_0 + n\,h , \quad n=0,1,2,\ldots . \) The fourth order Runge--Kutta method is based on computing yn+1 as follows
I. Classical RK4
By far the most often used is the classical fourth-order Runge-Kutta formula, which has a certain sleekness of organization about it:
The Runge-Kutta method iterates the x-values by simply adding a fixed step-size of h at each iteration. The y-iteration formula is far more interesting. It is a weighted average of four coefficients. Doing thi,s we see that k1 and k4 are given a weight of 1/6 in the weighted average, whereas k2 and k3 are weighted 1/3, or twice as heavily as k1 and k4. (As usual with a weighted average, the sum of the weights 1/6, 1/3, 1/3 and 1/6 is 1.) So what are these ki values that are being used in the weighted average?
- k1 is simply Euler's prediction.
- k2 is evaluated halfway across the prediction interval and gives us an estimate of the slope of the solution curve at this halfway point.
- k3 has a formula which is quite similar to that of k2 except that where k1 used to be, there is now a k2. Essentially, the f-value here is yet another estimate of the slope of the solution at the "midpoint" of the prediction interval. This time, however, the y-value of the midpoint is not based on Euler's prediction, but on the y-jump predicted already with k2.
- k4 evaluates f at \( x_n + h , \) which is at the extreme right of the prediction interval. The y-value coupled with this \( y_n + k_3 , \) is an estimate of the y-value at the right end of the interval, based on the y-jump just predicted by k3.
function R=rk4(f,a,b,ya,M)
% Input: f is the slope function entered as a string 'f'
% a and b are the left and right end points
% ya is the initial condition y(a)
% M is the number of steps
% Output: R=[T',Y'] where T is the vector of abscissas and Y is the vector of ordinates
h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=a:h:b;
Y(1)=ya;
for j=1:M
k1=h*feval(f,T(j),Y(j));
k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
k3=h*feval(f,T(j),Y(j));
k4=h*feval(f,T(j)+h,Y(j)+k3);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end
R=[T',Y'];
function rungekutta
h = 0.5; t = 0; w = 0.5;
fprintf(’Step 0: t = %12.8f, w = %12.8f\n’, t, w);
for i=1:4
k1 = h*f(t,w);
k2 = h*f(t+h/2, w+k1/2);
k3 = h*f(t+h/2, w+k2/2);
k4 = h*f(t+h, w+k3);
w = w + (k1+2*k2+2*k3+k4)/6;
t = t + h;
fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w);
end
%%%%%%%%%%%%%%%%%%
function v = f(t,y)
v = y^2-tˆ2+1;
II. Kutta Algorithm RK4
Kutta's algorithm (1901) of order four:
function Complete
III. Gill's Algorithm RK4
Gill's method:
function Complete
IV. Fifth order Butcher's method
Here is the fifth order Runge--Kutta method proposed by Butcher:
function Complete
V. Runge--Kutta--Fehlberg Method (RKF45)
One way to guarantee accuracy in the solution of an initial value problem is t solve the problem twice using step sizes h and h/2 and compare answers at the mesh points corresponding to the largest step size. However, this requires a significant amount of computation for the smaller step length and must be repeated if it is determined that the agreement is not good enough.
The Runge--Kutta--Fehlberg method (denoted RKF45) or Fehlberg method was developed by the German mathematician Erwin Fehlberg (1911--1990) in 1969 NASA report. The novelty of Fehlberg's method is that it is an embedded method from the Runge-Kutta family, and it has a procedure to determine if the proper step size h is being used. At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step length is reduced. If the answers agree to more significant digits than required, the step size is increased. RKF45 method is a method of order \( O( h^4 ) \) with an error estimator of order \( O( h^5 ) . \)Initially, Erwin Fehlberg introduced a so-called six stage Runge-Kutta method that requires six function evaluations per step. These function values can be combined with one set of coefficients to produce a fifth-order accurate approximation and with another set of coefficients to produce an independent fourth-order accurate approximation. Comparing these two approximations provides an error estimate and resulting step size control. Notice that it takes six stages to get fifth order. It is not possible to get fifth order with only five function evaluations per step. The new ode45 introduced in the late 1990s is based on an algorithm of Dormand and Prince (DOPRI). It uses six stages, employs the FSAL strategy, provides fourth and fifth order formulas, has local extrapolation and a companion interpolant. The new ode45 is so effective that the interpolant is called up to produce more finely spaced output at no additional cost measured in function evaluations. The codes for the two routines ode23 and ode45 are very similar.
Butcher tableau for Fehlberg's 4(5) method
Other Runge--Kutta methods could be found on
https://en.wikipedia.org/wiki/List_of_Runge-Kutta_methods
Some other popular 4/5 versions of Runge--Kutta methods:
- The Dormand--Prince 4/5 pair (Dormand, J. R.; Prince, P. J. (1980), "A family of embedded Runge--Kutta formulae", Journal of Computational and Applied Mathematics, 6 (1): 19–26, doi:10.1016/0771-050X(80)90013-3;
Dormand, John R. (1996), Numerical Methods for Differential Equations: A Computational Approach, Boca Raton: CRC Press, pp. 82--84, ISBN 0-8493-9433-3.). - The Cash--Karp pair. It was proposed by Professor Jeff R. Cash from Imperial College London and Alan H. Karp from IBM Scientific Center. (J. R. Cash, A. H. Karp. "A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides", ACM Transactions on Mathematical Software 16: 201--222, 1990. doi:10.1145/79505.79507).
- The Tsitorous 4/5 (http://users.ntua.gr/tsitoura/RK54_new_v2.pdf).
function R=rk45(f,a,b,ya,M,tol)
% Input: f is the slope function entered as a string 'f'
% a and b are the left and right end points
% ya is the initial condition y(a)
% M is the number of steps
& tol is the tolerance
% Output: R=[T',Y'] where T is the vector of abscissas and Y is the vector of ordinates
h=(b-a)/M;
big=1e15;
hmin=h/64;
hmax=64*h;
Y(1)=ya;
T(1)=a;
j=1;
br=b-0.00001*abs(b);
% Enter coefficients:
a2=1/4; b2=1/4; a3=3/8; b3=3/32; c3=9/32; a4=12/13;
b4=1932/2197; c4=-7200/2197; d4=7296/2197; a5=1;
b5=439/216; c5=-8; d5=3680/513; e5=-845/4104; a6=1/2;
b6=-8/27; c6=2; d6=-3544/2565; e6=1859/4104;
f6=-11/40; r1=1/360; r3=-128/4275; r4=-2197/75240; r5=1/50;
r6=2/55; n1=25/216; n3=1408/2565; n4=2197/4104; n5=-1/5;
while(T(j),b)
if ((T(j)+h) > br)
h=b-T(j);
end
k1=h*feval(f,T(j),Y(j));
y2=Y(j)+b2*k1;
if big < abs(y2)break,end
k2=Y(j)+b3*k1+c3*k2;
y3=Y(j)+b3*k1+c3*k2;
if big < abs(y3)break,end
k3=h*feval(f,T(j)+a3*h,y3);
y4=Y(j)+b4*k1+c4*k2+d4*k3;
if big < abs(y4)break,end
k4=h*feval(f,Y(j)+a4*h,y4);
y5=Y(j)+b5*k1+c5*k2+d5*k3+e5*k4;
if big < abs(y5)break,end
k6=h*feval(f,Y(j)+a6*h,y6);
err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
ynew=Y(j)+n1*k1+n3*k3+n4*k4+n5*k5;
% error and step size control
if((err < tol)|(h < 2*hmin))
Y(j+1)=ynew;
if((T(j)+h) > br)
T(j+1)=b;
else
T(j+1)=T(j)+h;
end
j=j+1;
end
if (err==0)
s=0;
else
s=0.84*(tol*h/err)^(0.25);
end
if((s < 0.75)&(h > 2*hmin))
h=h/2;
end
if((s > 1.50)&(2*h < hmax))
h=2*h
end
if((big < abs(Y(j)))|(max1==j)),break,end
M=j;
if (b > T(j))
M=j+1;
else
M=j;
end
end
R=[T',Y'];
function rk45
epsilon = 0.00001;
h = 0.2; t = 0;
w = 0.5; % initial value
i = 0;
fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w);
while t < 2
h = min(h, 2-t);
k1 = h*f(t,w);
k2 = h*f(t+h/4, w+k1/4);
k3 = h*f(t+3*h/8, w+3*k1/32+9*k2/32);
k4 = h*f(t+12*h/13, w+1932*k1/2197-7200*k2/2197+7296*k3/2197);
k5 = h*f(t+h, w+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
k6 = h*f(t+h/2, w-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
w1 = w + 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
w2 = w + 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55;
R = abs(w1-w2)/h;
delta = 0.84*(epsilon/R)ˆ(1/4);
if R <= epsilon
t = t+h;
w = w1;
i = i+1;
fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w);
h = delta*h;
else h = delta*h;
end
end
%%%%%%%%%%%%%%%%%%
function v = f(t,y)
v = y^2-tˆ2+1;
function RungeKuttaMethod
%example of Runge-Kutta method of order 4.
%y' = y - t^2 + 1.
a = 0; %time interval we are solving over
b = 4;
N = 1000; %number of steps
t = zeros(1, N); %preallocate memory for these variables - this is preferred but is not required to work%
w = zeros(1, N);
w(1) = 0.5; %initial value of w.
h = (b - a)/N; %step size
t(1) = a; %initial vlaue for time.
F = @(t, y) y - t^2 + 1; %thing we are solving.
for i = 1:(N-1)
K1 = h*(F(t(i), w(i))); %this is the RungKutta Method.
K2 = h*(F(t(i) + 0.5*h, w(i) + 0.5*K1));
K3 = h*(F(t(i) + 0.5*h, w(i) + 0.5*K2));
K4 = h*(F(t(i) + h, w(i) + K3));
w(i+1) = w(i) + (K1 + 2*K2 + 2*K3 + K4)/6;
t(i+1) = a + i*h;
end
%Fa = (t+1).^2 -0.5*exp(t); %actualy solution we ca use to compare between both methods%
plot(t, w)
%hold on %use these 3 lines to see both solutions plotted on
%plot(t, Fa, 'r') %same graph for comparison.
%hold off
end