Page 278 - Applied Numerical Methods Using MATLAB

RUNGE–KUTTA METHOD 267

function [t,y] = ode_Heun(f,tspan,y0,N)
%Heun method to solve vector differential equation y’(t) = f(t,y(t))
% for tspan = [t0,tf] and with the initial value y0 and N time steps
if nargin<4|N<=0, N= 100; end
if nargin<3, y0 = 0; end
h = (tspan(2) - tspan(1))/N; %stepsize
t = tspan(1)+[0:N]’*h; %time vector
y(1,:) = y0(:)’; %always make the initial value a row vector
for k = 1:N
fk = feval(f,t(k),y(k,:)); y(k+1,:) = y(k,:)+h*fk; %Eq.(6.2.3)
y(k+1,:) = y(k,:) +h/2*(fk +feval(f,t(k+1),y(k+1,:))); %Eq.(6.2.4)
end

But, the right-hand side (RHS) of this equation has y k+1 , which is unknown at
t k . To resolve this problem, we replace the y k+1 on the RHS by the following
approximation:
∼
y k+1 = y k + hf(t k , y k ) (6.2.3)
so that it becomes

h
y k+1 = y k + {f(t k , y k ) + f(t k+1 , y k + hf(t k , y k ))} (6.2.4)
2
This is Heun’s method, which is implemented in the MATLAB routine
“ode_Heun()”. It is a kind of predictor-and-corrector method in that it predicts
the value of y k+1 by Eq. (6.2.3) at t k and then corrects the predicted value by
2
Eq. (6.2.4) at t k+1 . The truncation error of Heun’s method is O(h ) (proportional
2
to h ) as shown in Eq. (5.6.1), while the error of Euler’s method is O(h).

6.3 RUNGE–KUTTA METHOD

Although Heun’s method is a little better than the Euler’s method, it is still not
accurate enough for most real-world problems. The fourth-order Runge–Kutta
4
(RK4) method having a truncation error of O(h ) is one of the most widely used
methods for solving differential equations. Its algorithm is described below.
h
y k+1 = y k + (f k1 + 2f k2 + 2f k3 + f k4 ) (6.3.1)
6
where

f k1 = f(t k , y k ) (6.3.2a)
f k2 = f(t k + h/2, y k + f k1 h/2) (6.3.2b)
f k3 = f(t k + h/2, y k + f k2 h/2) (6.3.2c)
f k4 = f(t k + h, y k + f k3 h) (6.3.2d)