Table Of ContentNumerical Methods for Differential Equations
YA YAN LU
Department of Mathematics
City University of Hong Kong
Kowloon, Hong Kong
1
Contents
1 ODEIVP:ExplicitOne-stepMethods 4
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 EulerandRunge-KuttaMethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Localtruncation errorandorder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4 EmbeddedRunge-Kuttamethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 ODEIVP:ImplicitOne-stepMethods 17
2.1 Stiffequations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Implicitone-step methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3 ODEIVP:Multi-stepMethods 23
3.1 Explicitmulti-step methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 Implicitmulti-step methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 ODEIVP:StabilityConcepts 29
4.1 Zerostability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Absolutestability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5 ODEBoundaryValueProblems 34
5.1 Theshootingmethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.2 Finitedifference methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3 Thefiniteelementmethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6 FiniteDifferenceMethodsforParabolicPDEs 42
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.2 Classicalexplicit method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6.3 Crank-Nicolson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
6.4 Stabilityanalysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.5 Alternating direction implicitmethod . . . . . . . . . . . . . . . . . . . . . . . . . . 48
7 FiniteDifferenceMethodsforHyperbolicPDEs 52
7.1 Firstorderhyperbolic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
7.2 Explicitmethodsforwaveequation . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.3 Maxwell’sequations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2
8 FiniteDifferenceMethodsforEllipticPDEs 60
8.1 Finitedifference methodforPoissonequation . . . . . . . . . . . . . . . . . . . . . . 60
8.2 FastPoissonsolverbasedonFFT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8.3 Classicaliterativemethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
8.4 Conjugategradient method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.4.1 1-Doptimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.4.2 Subspaceminimization problem . . . . . . . . . . . . . . . . . . . . . . . . . 65
8.4.3 Orthogonal residual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
8.4.4 Thenextconjugate direction . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
8.4.5 Themethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
8.4.6 Rateofconvergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3
Chapter 1
ODE IVP: Explicit One-step Methods
1.1 Introduction
In this chapter, we study numerical methods for initial value problems (IVP) of ordinary differential
equations (ODE).Thefirststepistore-formulate yourODEasasystemoffirstorderODEs:
dy
= f(t,y) for t >t (1.1)
0
dt
withtheinitialcondition
y(t )=y (1.2)
0 0
wheretistheindependentvariable,y=y(t)istheunknownfunctionoft,y isthegiveninitialcondition,
0
and f is a given function of t and y which describes the differential equation. High order differential
equationscanalsobewrittenasafirstordersystembyintroducingthederivativesasnewfunctions. Our
numerical methods can be used to solve any ordinary differential equations. We only need to specify
thefunction f.
Thevariablet isdiscretized, sayt for j=0,1,2,..., thenwedeterminey y(t )for j=1,2,3,....
j j j
≈
Thefirstclassofmethods(Runge-Kuttamethods)involveone-step. Ify iscalculated,thenweconstruct
j
y fromy . Previousvaluessuchasy arenotneeded. SincethisisanIVPandforthefirststep,we
j+1 j j 1
−
have y only att , then wecan findy , y , ...,inasequence. Theone-step methods are varynatural. A
0 0 1 2
higherordermethodgivesamoreaccuratenumericalsolutionthanalowerordermethodforafixedstep
size. Butahigher orderone-step methodrequires moreevaluations ofthe f function. Forexample, the
firstorderEuler’smethodrequiresonlyoneevaluationof f,i.e., f(t ,y ),butafourthorderRunge-Kutta
j j
methodrequires fourevaluations of f.
For a large scale problem, the computation of f could be time consuming. Thus, it is desirable to
havehighordermethodsthatrequireonlyoneevaluationof f ineachstep. Thisisnotpossibleinaone-
stepmethod. Butitispossible inamulti-step method. Therefore, themainadvantage ofthemulti-step
methodisthattheyareefficient. However,theyaremoredifficulttouse.
Forone-step methods, wewillintroduce implicit methods. These aremethods designed forthe so-
called“stiff”ODEs. Ifanexplicit methodisusedforastiffODEandthestepsizeisnotsmallenough,
the error (between the exact and the numerical solution) maygrow very fast. Forthese stiff ODEs,the
4
implicit methods are useful. The situation is the same for multi-step methods. We also need implicit
multi-stepmethodsforstiffODEs.
We will also introduce the embedded Runge-Kutta methods. These are methods that combine two
methods together, so that the step size can be automatically chosen for a desired accuracy. There are
also multi-step methods that allow automatic selection of the step size. Butthey are more complicated
andwewillnotcoverthem.
Considerthefollowingexample. Wehavethefollowingdifferential equation foru=u(t):
u
u +sin(t) 1+(u )2u + =t2 (1.3)
′′′ q ′′ ′ 1+e−t
fort >0,withtheinitialconditions:
u(0)=1, u(0)=2, u (0)=3. (1.4)
′ ′′
Wecanintroduce avectory
u(t)
y(t)=u′(t)
u (t)
′′
andwritedowntheequation foryas
u
′
y′= f(t,y)= u′′
sin(t) 1+(u )2 u u/(1+e t)+t2
′′ ′ −
− −
p
Theinitialcondition isy(0)=[1,2,3]. HereisasimpleMATLABprogramfortheabovefunction f.
function k = f(t, y)
% remember y is a column vector of three components.
k = zeros(3,1);
k(1) = y(2);
k(2) = y(3);
k(3) = -sin(t) * sqrt(1+y(3)^2) * y(2) - y(1)/(1 + exp(-t)) + t^2;
In the MATLAB program, y(1), y(2), y(3) are the three components of the vector y. They are
u(t), u(t) and u (t), respectively. They are different from y(1), y(2) and y(3) which are the vectors y
′ ′′
evaluatedatt =1,t =2andt =3. Noticethatwealsohavey(0),whichistheinitialvalueofy. Butwe
donothavey(0). Anyway,thecomponents ofyareonlyusedinsidetheMATLABprograms.
A numerical method is usually given for the general system (1.1-1.2). We specify the system of
ODEsbywriting aprogram for thefunction f, then the samenumerical method can beeasily used for
solvingmanydifferentdifferential equations.
1.2 Euler and Runge-Kutta Methods
Numericalmethodsstartwithadiscretization oft byt ,t ,t ,...,say
0 1 2
t =t + jh
j 0
5
wherehisthestepsize. Numericalmethodsareformulasfory ,y ,y ,...,wherey istheapproximate
1 2 3 j
solution att . Weusey(t )todenotethe(unknown)exactsolution, thus
j j
y y(t ).
j j
≈
Pleasenoticethatwhenyisavector,y ,y ,...,arealsovectors. Inparticular,y isnotthefirstcomponent
1 2 1
ofyvector, y isnot the 2nd component of theyvector. Thecomponents ofyare only explicitly given
2
insidetheMATLABprogramsasy(1), y(2),etc.
Euler’smethod:
y =y +hf(t ,y ). (1.5)
j+1 j j j
Sincey istheknowninitialcondition,theaboveformulaallowsusofindy ,y ,etc,inasequence. The
0 1 2
Euler’s method can be easily derived as follows. First, we assume h is small and consider the Taylor
expansion:
y(t )=y(t +h)=y(t )+hy(t )+...
1 0 0 ′ 0
Now, we know that y(t )= f(t ,y(t )). If we keep only the first two terms of the Taylor series, we
′ 0 0 0
obtainthefirststepofEuler’smethod:
y =y +hf(t ,y ),
1 0 0 0
wherey(t )isreplaced bythe“numerical solution” y ,etc. Thegeneralstepfromt tot issimilar.
1 1 j j+1
HereisaMATLABprogram fortheEuler’smethod:
function y1 = eulerstep(h, t0, y0)
% This is one step of the Euler’s method. It is
% given for the first step, but any other step
% is just the same. You need the MATLAB function
% f to specify the system of ODEs.
y1 = y0 + h* f(t0, y0)
Now,letussolve(1.3-1.4)fromt=0tot=1withthestepsizeh=0.01. Forthispurpose,weneed
towriteamainprogram. Inthemainprogram,wespecifytheinitialconditions,initialtimet ,finaltime
0
andthetotalnumberofsteps. Thestepsizecanthenbecalculated. HereistheMATLABprogram.
% The main program to solve (1.3)-(1.4) from t=0 to
% t = 1 by Euler’s method.
% initial time
t0 = 0;
% final time
tfinal = 1;
% number of steps
nsteps = 100;
% step size
6
h = (tfinal - t0)/ nsteps;
% initial conditions
y = [1, 2, 3]’;
% set the variable t.
t = t0
% go through the steps.
for j= 1 : nsteps
y = eulerstep(h, t, y)
t = t + h
% saved output for u(t) only, i.e. the first component of y.
tout(j) = t;
u(j) = y(1);
end
% draw a figure for the solution u.
plot(tout, u)
Now,insiderMATLAB,inafoldercontainingthethreeprograms: f.m, eulerstep.m, eulermain.m,
ifwetype eulermain,wewill see asolution curve. That isthe solid curve in Fig. 1.1. Thisis for the
Numerical solutions by Euler’s method using h=0.01, 0.1, 0.2
5
4.5
4
3.5
u(t) 3
2.5
2
1.5
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
Figure 1.1: Numerical solutions of (1.3) and (1.4) by Euler’s method. The solid curve is for h=0.01.
The“+”isforh=0.1andthe“o”isforh=0.2.
case of h=0.01. We also want to see what happens if h is 0.2 and 0.1. For this purpose, we change
nsteps to 5 and 10, then use plot(tout, u, ’o’) and plot(tout, u, ’+’) to show the results.
AllthreeplotsareshownintheFig. 1.1.
TheEuler’smethodisnotveryaccurate. Toobtainanumericalsolutionwithanacceptableaccuracy,
wehavetouseaverysmallstepsizeh. Asmallstepsizehimpliesalargernumberofsteps, thusmore
7
computing time. It is desirable to develop methods that are more accurate than Euler’s method. If we
lookattheTaylorseriesagain,wehave
h2 h3
y(t )=y(t +h)=y(t )+hy(t )+ y (t )+ y (t )+...
1 0 0 ′ 0 ′′ 0 ′′′ 0
2 6
Thiscanbewrittenas
y(t ) y(t ) h h2
1 0
− =y′(t0)+ y′′(t0)+ y′′′(t0)+... (1.6)
h 2 6
Actually, therighthandsideisamoreaccurateapproximation fory(t +h/2),since
′ 0
h h h2
y(t + )=y(t )+ y (t )+ y (t )+...
′ 0 ′ 0 ′′ 0 ′′′ 0
2 2 8
The first two terms on the right hand sides of the above two equations are identical, although the third
termsinvolving y (t )aredifferent. Thus,
′′′ 0
y(t ) y(t ) h h h
1 0
− y′(t0+ )= f t0+ ,y(t0+ )
h ≈ 2 (cid:18) 2 2 (cid:19)
The right hand side now involves y(t +h/2). Of course, this is now known, because we only have
0
y(t ). The idea is that we can use Euler’s method (with half step size h/2) to get an approximate
0
y(t +h/2),thenusetheabovetogetanapproximationofy(t ). TheEulerapproximationfory(t +h/2)
0 1 0
isy(t )+h/2f(t ,y ). Therefore, wehave
0 0 0
k = f(t ,y ) (1.7)
1 0 0
h h
k = f(t + ,y + k ) (1.8)
2 0 0 1
2 2
y = y +hk . (1.9)
1 0 2
Thisisthefirststepoftheso-calledmidpointmethod. Thegeneralstepisobtainedbysimplyreplacing
t ,y andy byt ,y andy ,respectively.
0 0 1 j j j+1
Therighthandsideof(1.6)canalsobeapproximated by(y(t )+y(t ))/2,because
′ 0 ′ 1
y(t )+y(t ) h h2
′ 0 ′ 1
=y(t )+ y (t )+ y (t )+...
′ 0 ′′ 0 ′′′ 0
2 2 4
Therefore, wehave
y(t ) y(t ) y(t )+y(t )
1 0 ′ 0 ′ 1
− .
h ≈ 2
We can replace y(t ) and y(t ) by f(t ,y(t )) and f(t ,y(t )), but of course, we do not know y(t ),
′ 0 ′ 1 0 0 1 1 1
becausethatiswhatwearetryingtosolve. ButwecanuseEuler’smethodtogetthefirstapproximation
ofy(t )anduseitin f(t ,y(t )),thenusetheabovetogetthesecond(andbetter)approximationofy(t ).
1 1 1 1
Thiscanbesummarizedas
k = f(t ,y ) (1.10)
1 0 0
k = f(t +h,y +hk ) (1.11)
2 0 0 1
h
y = y + (k +k ). (1.12)
1 0 1 2
2
8
Thisisthefirststepoftheso-calledmodifiedEuler’smethod. Thegeneralstepfromt tot iseasily
j j+1
obtained byreplacing thesubscripts 0and1by jand j+1,respectively.
Similarly,therighthandsideof(1.6)canbeapproximated by
Ay(t )+By(t +a h),
′ 0 ′ 0
wherea isagivenconstant, 0<a 1,thecoefficients AandBcanbedetermined, suchthattheabove
≤
matchesthefirsttwotermsoftherighthandsideof(1.6). Weobtain
1 1
A=1 , B= .
−2a 2a
Theny(t +a h)= f(t +a h,y(t +a h))andweuseEuler’smethodtoapproximate y(t +a h). Thatis
′ 0 0 0 0
y(t +a h) y(t )+a hf(t ,y(t )).
0 0 0 0
≈
Finally,weobtainthefollowinggeneral2ndorderRunge-KuttaMethods:
k = f(t ,y ) (1.13)
1 j j
k = f(t +a h,y +a hk ) (1.14)
2 j j 1
1 1
y = y +h 1 k + k (1.15)
j+1 j (cid:20)(cid:18) −2a (cid:19) 1 2a 2(cid:21)
Since a is an arbitrary parameter, there are infinitely many2nd order Runge-Kutta methods. Themid-
point method and the modified Euler’s method correspond to a =1/2 and a =1, respectively. In this
formula,k andk aretemporaryvariables, theyaredifferent fordifferent steps.
1 2
There aremanyother Runge-Kutta methods (3rd order, 4thorder andhigher order). Thefollowing
classical 4thorderRunge-Kuttamethodiswidelyused,because itisquiteeasytoremember.
k = f(t ,y ) (1.16)
1 j j
h h
k = f(t + ,y + k ) (1.17)
2 j j 1
2 2
h h
k = f(t + ,y + k ) (1.18)
3 j j 2
2 2
k = f(t +h,y +hk ) (1.19)
4 j j 3
h
y = y + (k +2k +2k +k ) (1.20)
j+1 j 1 2 3 4
6
Wehavementionedtheorderofamethodabove. Thisconceptwillbeexplained inthenextsection.
Next,weconsider aMATLABimplementation ofthemidpoint method. Forthispurpose, wewrite
thefollowingfunction calledmidptstepwhichissavedinthefilecalledmidptstep.m.
function y1 = midptstep(h, t0, y0)
% This is midpoint method (one of the second order Runge-Kutta methods).
% It is given for the first step, but any other step is just the same.
% You need the MATLAB function f to specify the system of ODEs.
k1 = f(t0, y0);
k2 = f(t0+h/2, y0 + (h/2)*k1)
y1 = y0 + h* k2;
9
To solve the same differential equation (1.3-1.4), we need the earlier MATLAB function f and a
main program. We can write a main program by copying the main program eulermain for Euler’s
method. The new mainprogram midptmainisdifferent from eulermainonly in one line. Theline y
= eulerstep(h, t, y)isnowreplacedby
y = midptstep(h, t, y)
You can see that writing a program for a new method is very easy, since we have separated the differ-
ential equation (inf.m)andthenumerical method(ineulerstep.mormidptstep.m)fromthemain
program. In Fig. 1.2, we show the numerical solution u(t) for (1.3-1.4) calculated by the midpoint
Numerical solutions by the midpoint method for h=0.01 and h=0.2
5
4.5
4
3.5
u(t) 3
2.5
2
1.5
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
Figure1.2: Numericalsolutionsbythemidpointmethod. Thesolidcurveisforh=0.01. The“o”isfor
h=0.2.
method with h= 0.01 and h=0.2. You can see that the midpoint solution obtained with h=0.2 is
muchmoreaccurate thantheEuler’ssolution withthesameh.
1.3 Local truncation error and order
Whenanumericalmethodisusedtosolveadifferential equation, wewanttoknowhowaccurate isthe
numerical solution. Wewilldenote the exact solution as y(t), thus y(t )is the exact solution att . The
j j
numericalsolution att isdenotedasy ,therefore, weareinterested inthefollowingerror:
j j
e = y(t ) y .
j j j
| − |
We do not expect to be able to know e exactly, because we do not have the exact solution in general.
j
Therefore, wewillbehappy tohave someestimates (such asapproximate formulas orinequalities) for
e . However,eventhisisnotsoeasy. Thereason isthattheerroraccumulates. Letuslookatthesteps.
j
10
