COMPACT NUMERICAL METHODS FOR STIFF
DIFFERENTIAL EQUATIONS
By
EDWARD PURBA
Sarjana S-1
Bandung Institute ofTechnology
Bandung, Indonesia
1984
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
May, 1996
COMPACT NUMERICAL METHODS FOR STIFF
DIFFERENTIAL EQUATIONS
Thesis Approved:
8Th~
f-1. vJ
fJJtJ c /t 1,-i //7 ;/·
~' ~) a/ ·trPA. / . ?kdJ.: \L/ ·/
.:iAtJYntM (!. ~
Dean of the Graduate College
11
ACKNOWLEDGMENTS
First, I would like to express my sincere appreciation to my major adviser, Dr.
John P. Chandler, for providing continuous support in finishing this thesis, and helping
me with his intelligent supervision, constructive guidance, and insightful critiques and
suggestions from the early until the late stages of the works. My sincere appreciation is
also addressed to my other committee members Dr. Blayne E. Mayfield and Dr. Huizhu
Lu, for their cooperation, friendship, and assistance.
So many people have helped me in completing this thesis, it is impossible to
acknowledge them all personally. I have gained much value from my relationship with H.
Sinaga during my study at Oklahoma State University, and encouragement from my
friends including Y. Sofyan, T. Rahardjo and Dr. Robert A. Divali.
I also am indebted to the Government of Indonesia, who sponsored my graduate
study through the Ministry for Research and Technology, in the framework of the
implementation of the Science And Technology for Industrial Development. Many
thanks to Prof. Dr. Ing. B.J. Habibie, State Minister of Research and Technology, Prof. S.
Sapiie D.Sc., Mrs. Ir. Ina Juniarti, Vice President for IPTN Computing Center, and Mr.
Ir. Muriawan A. Kadir.
Finally, I would like to dedicate this thesis to my mother, late father and late
sister Ida Purba.
iii
TABLE OF CONTENTS
Chapter Page
I. INTitOIJlJCTION --------------------------------------------------------------------- 1
1.1 Initial Value Problem for First-Order OIJEs --------------------------------- 2
1.2 The Existence and lJniqueness of Solution of Initial Value Problems --- 3
1.3 Taylor's Series ------------------------------------------------------------------- 4
1.4 Stability and Instability of OIJEs ---------------------------------------------- 6
1.5 Convergence of Euler Method ------------------------------------------------- 9
1.6 The Multistep Concept ---------------------------------------------------------- 12
II. STABILITY ANALYSIS ------------------------------------------------------------- 16
2.1 Stability of Euler Methods ----------------------------------------------------- 18
2.2 Stability of Explicit Itunge-Kutta Methods ---------------------------------- 24
2.3 Stability oflmplicit Itunge-Kutta Methods ---------------------------------- 30
2.4 l)ifference Equations------------------------------------------------------------- 31
III. STABILITY OF MlJL TISTEP METHOIJS --------------------------------------- 35
3.1 Linear l)ifference Equations ---------------------------------------------------- 3 8
3 .2 Adams Methods ------------------------------------------------------------------ 44
3.3 Backward l)ifferentiation Formula (BIJF) ----------------------------------- 48
3 .4 Predictor-Corrector Methods --------------------------------------------------- 51
3.5 Extrapolation Methods for Solving OIJEs ----------------------------------- 53
IV. STIFF OWINAitY IJIFFEitENTIAL EQlJATIONS --------------------------- 62
4.1 Stiff l)ifferential Problems ----------------------------------------------------- 66
4.2 Stiffness Concepts --------------------------------------------------------------- 68
4.3 Stability Concepts ofNumerical Stiffness Methods ------------------------ 72
4.4 A Modified Euler Method for Solving OIJEs ------------------------------- 75
4.5 An Explicit Exponential Method for Solving OIJEs ----------------------- 79
4. 6 0 l) E Solvers ---------------------------------------------------------------------- 83
V. ANALYSIS ANI) IJISClJSSION --------------------------------------------------- 86
5.1 Implementation of a Modified Euler Method-------------------------------- 87
iv
Chapter Page
5.2 Implementation of an Exponential Method ---------------------------------- 90
5.3 Analysis of Kidney Problems -------------------------------------------------- 92
5.4 Analysis of Autocatalytic Problems ------------------------------------------ 95
5.5 Analysis of Problem D4 of Enright et al. ------------------------------------ 97
5.6 Analysis of Gupta and Wallace's Problem----------------------------------- 99
VI. CONCLUSION AND SUGGESTION --------------------------------------------- 101
SELECTED BIBLIOGRAPHY ------------------------------------------------------------ 103
APPENDIXES -------------------------------------------------------------------------------- 116
APPENDIX A-- SOFTWARE AVAILABILITY -------------------------------- 116
APPENDIX B --ROOTS OF A COMPLEX POLYNOMIAL------------------ 120
APPENDIX C-- LINEAR MULTISTEP FORMULAS ------------------------- 128
APPENDIX D-- ESTIMATION OF ERROR AND STEPSIZE
CONTROL --------------------------------------------------------- 13 7
APPENDIX E -- STABILITY REGION OF SIMPLE PREDICTORCORRECTOR
METHODS ------------------------------------- 148
APPENDIX F --COLLECTION OF PROGRAMS ------------------------------ 152
APPENDIX G-- COLLECTION OF TABLES ----------------------------------- 174
v
LIST OF TABLES
Table Page
2.1. Stability Region of Runge-Kutta for Orders 1 to 10 ------------------------------- 28
2.2 Butcher's Array ofRunge-Kutta Methods ------------------------------------------ 29
2.3 Butcher's Array of a Fourth-Order Runge-Kutta method------------------------- 30
3.1 A Trial Solution to Find a Particular Solution ofNonhomogeneous
Difference Equations ------------------------------------------------------------------- 42
3 .2 Table of Extrapolation------------------------------------------------------------------ 56
3.3 Romberg Representation ofy'= -y, y(O) = 1 ---------------------------------------- 59
3.4 The Computation Result ofy'= -y, y(0)=1, Using Extrapolation,
Euler, and ABM-4 Methods----------------------------------------------------------- 60
G.1 Table ofPerformances ofMEBDF, VODE, LSODE, and EPSODE
for the Case of Kidney Problem ------------------------------------------------------ 17 4
G.2 Table of Performances ofMEBDF, VODE, LSODE, and EPSODE
for the Case of Autocatalitic Reaction Pathway Problem------------------------- 184
G.3 Table of Performances ofMEBDF, VODE, LSODE, and EPSODE
for the Case of D4 of Enright et al. --------------------------------------------------- 186
G.4 Table of Performances ofMEBDF, VODE, LSODE, and EPSODE
for the Case of problem proposed by Gupta and Wallace------------------------- 188
G.5 The Computation Results for Problem 1 of Chapter 4.4
Using Modified Euler Method and Exact Solution -------------------------------- 190
G.6 The Computation Results for Problem 2 of Chapter 4.4
Using Modified Euler Method and Mathematica Packages----------------------- 191
vi
Table Page
G.7 The Computation Results for Problem 1 of Chapter 4.4
Using Exponential Method and the Exact Solution ------------------------------- 192
G.8 The Computation Results for Problem 2 of Chapter 4.4
Using Exponential Method and Mathematica Packages-------------------------- 193
vii
LIST OF FIGURES
Figure Page
1.1 Solution y = (3+ 3E/5)e-21 + 2E/5e31 forE= 0, 0.01, and 0.1----------------------- 7
1.2 Crude Euler Method -------------------------------------------------------------------- 10
1.3 Numerical Results Using Exact, Euler, and Trapezoidal Methods-------------- 11
2.1 Stability Region of the Explicit Euler Method ------------------------------------- 19
2.2 Stability Region of the Implicit Euler Method ------------------------------------- 19
2.3 Euler Method with Stepsize h = 0.01, 0.03, 0.04, and 0.05 ---------------------- 21
2.4 Crude Euler, Implicit Euler Using PC, and Implicit Euler Using Newton's --- 23
2.5 Implicit Euler by Showing the Vector Gradient ----------------------------------- 24
2. 6 A Second Order of Runge-Kutta ----------------------------------------------------- 26
2.7 Stability Region of Second order Runge-Kutta ------------------------------------ 27
2.8 Difference Equation Methods Using Forward, Backward and Central
for h = 0.05, 1, 2, and 3.---------------------------------------------------------------- 34
3.1 Computational Results ofyn+2 = -4Yn+J+ 5yn+ h(4fn+J+ 2fn)
Applied toy' = y, y(O) = 1------------------------------------------------------------- 43
3.2 Computation Results ofy' = -y, y(O) = 1 Using Extrapolation,
Euler, and ABM -4 Methods----------------------------------------------------------- 61
4.1 Graphical Solutions of e-t and e-Joot -------------------------------------------------- 64
4.2 A Representation of a Minimum Region of A-Stability -------------------------- 73
4.3 A Representation of a Minimum Region of A(a.)-Stability----------------------- 74
viii
Figure Page
4.4 A Representation of a Minimum Region of A0-Stability ------------------------ 74
4.5 A Representation of a Minimum Region of a Stiffly-Stable Method ---------- 75
5.1 Problem 1 of Section 4.4 Using Modified Euler Method ------------------------ 88
5.2 Problem 2 of Section 4.4 Using Modified Euler Method ------------------------ 89
5.3 Problem 1 of Section 4.4 Using Exponential Method ---------------------------- 90
5.4 Problem 2 of Section 4.4 Using Exponential Method---------------------------- 91
5.5 Kidney Problems with A= 902688359, 0.0990283499,
0. 992 5 2113 41 , and 1. 0 3 048 79 8 56 -------------------------------------------------- 92
5.6 Kidney Problems with A= 0.99, 0.9, and 0----------------------------------------- 93
5.7 Output ofMathematica for Robertson Problem------------------------------------ 95
5.8 Robertson Problem Solved with EPSODE ----------------------------------------- 96
5.9 Problem D4 of Enright et al. Solved with Mathematica ------------------------- 97
5.10 Problem D4 of Enright et al. Solved with MEBDF ------------------------------- 98
5.11 Gupta and Wallace's Problem Solved with Mathematica ----------------------- 99
5.12 Gupta and Wallace's Problem Solved with VODE ------------------------------- 100
B-1 Stability Region ofRunge-Kutta Methods Orders 1, 2, 3, 4, and 5------------- 127
E-1 A PC-Method Where a Trapezoidal as a Corrector and
an Euler as a Predictor ---------------------------------------------------------------- 150
E-2 A PC-Method Where the Implicit Euler Method as Corrector
and Crude Euler as a Predictor ------------------------------------------------------- 151
ix
LIST OF SYMBOLS AND TERMS
Symbols Meanings
< less
~ less or equal
> greater
~ greater or equal
* not equal
I I absolute value
such that
{ } set
c subset
00 infinity
mm minimum
max maximum
factorial
I summation
lim limit
x~b x approaches b
V' backward difference
X
Symbols
Re(z)
Im(z)
arg(z)
u
>>
a !at
0()
[ ] T
~
Meanings
the real part of a complex number z
the imaginary part of a complex number z
argument of a complex number z
union
far greater
partial derivative over t
order
matrix transpose
set of real numbers
xi
CHAPTER I
INTRODUCTION
Before computers were involved in human lives, scientific practice only
recognized two terminologies; theory and experiment. However, after computers become
involved in human lives, computational science becomes one of the terminologies in
scientific practice and stands beside of the others as an essential methodology. It is
undeniable that most parts of science and engineering can be modeled by mathematical
equations. Unfortunately, most of mathematical equations are not easy to solve exactly,
so we need to approximate their solutions. The methods used to approximate these
solutions are called numerical methods. Scientific computing starts by transforming the
real-life problems into mathematical models. In this case, essential features of scientific
problems are represented in the form of mathematical equations. The second step is to
modify the mathematical models, so that they are suitable for numerical computations. In
this stage, usually we have to discretizise the domain of problems into steps of
computation, so that the formulations can be executed numerically either by computers or
by human brains manually, using the operators +,-, *, and I. The third stage is to test and
validate the solutions. This can be done by: comparing with the exact solutions,
comparing with previously validated computations, comparing with laboratory
experiments, and analyzing the convergences, the errors, and the diagnostic data. Finally,
1
2
the tested and validated codes are ready to use for various scientific and engineering
problems.
Some scientific and engineering problems are modeled by ordinary differential
equations such as y'= f(t,y), y(t0) = y0. Among them are chemical reactions (chemistry),
heat-flow problems (thermodynamics), electrical circuits (electrical engineering), force
problems (mechanics), rate of bacterial growth (biological science), decompositions of
radioactive material (atomic physics), population growth (statistics), and simulation and
control systems (control engineering). Due to difficulties in finding the exact solutions of
ordinary differential equations (ODEs), the study of numerical solution of ODEs becomes
important.
The choice of numerical methods used for the approximations depends on the
accuracy needed and the characteristic of the problem. Since every nth-order of ODE can
be transformed into a system of first-order ODEs, here the discussion will be confined
only to systems of first-order ODEs. For a similar reason, the discussion will also be
limited to initial value problems (IPVs).
1.1 Initial Value Problem for First-Order ODEs
The general form of an initial value problem for a first-order ODE can be written
as y'= f(t,y), where an initial value y(a) is given and tis in the interval [a,b]. Every nthorder
ODE in the form of d"y/dt" = y<n) = f(t,y,y< 1>,y<2>, ... ,y<n-l)) with initial conditions
y(i)(a)= ci, i=0,1, ... ,n-1, can be converted into a system of first-order ODEs ofthe form
dy/dt = fi(t,y~>y2, ... ,yn), yiCa)=ci , i=1,2, ... ,n-1 by substituting y1=y, y2=dy/dt, ..
3
d 3Y (d y )2 1
.,y0=dy0 _/dt. For example, a third-order IPV - 3 + - + 3y = e , y(O) = 1, y'(O) = 0,
dt dt
y"(O) = 0 can be converted into an initial-value problem for the variables y, dy/dt, and
d2y/dt2 by substituting y1=y, y2=dy/dt, and y3=d2y/dt2
. Applying these new variables will
transform the third-order ODE into a system of first-order ODEs that can be written as
dy I = dy 2 - dy 3 - I 2
dt y2, dt -y3, dt -e -y2 -3yi
where y1(0)=0, y2(0)=0, yiO)=O. After transforming the problem into a system of first-order
ODEs, our problem becomes how to solve this system of first-order ODEs using
numerical computing methods. Solving the first-order ODE will lead us to nonstiff or
stiff differential systems. These problems affect the stability of the methods used to
solve the system of first-order ODEs.
1.2 The Existence and Uniqueness of Solutions of Initial Value Problems
The first interesting question in handling first-order ODEs is to ask about the
existence of a solution. By investigating the existence of a solution, we do not waste our
time trying to solve a problem that has no solution in a given domain. An understanding
of Lipschitz conditions plays an important rule in proving the existence and the
uniqueness of solutions of first-order ODEs.
A function f(t,y) is said to satisfy a Lipschitz condition with respect to y in a
region D c 9tx9t, if for all (t,y) and (t,z) in D, there exists a constant L > 0 such that
lf(t,y) - f(t,z)l ~ Lly-zl. Here, the constant L is called a Lipschitz constant for f. For
example, let f(t,y) = 1 + / for D={(t,y) I ltl < 1, IYI < 1}. It is easy to see that
4
2 2 2 2 . lf(t,y)-f(t,z)l = l1+y - (l+z )I= IY -z I= ly+zlly-zl. Smce ltl < 1 and IYI < 1, than we can
take L=2 as a Lipschitz constant for f.
Many authors in the field of Differential Equations have proved the existence of
solutions of first-order ODEs. Usually, a constructive solution of a first-order ODE can
be accomplished using Picard iteration. Braun [10] said that if f and 8f /8y are
continuous in the rectangle D={(t,y) I t0 ::; t::; t0+a, IY- Yol ::; b }, then the initial value
problem y'=f(t,y), y(t0)=y0 has at least one solution y(t) on the interval t0 ::; t ::; t0+a,
where a= min{a,b/M} and M = max lf(t,y)l. Shampine and Gordon [99] used the
(t,y)inD
terminology of Lipschitz condition to show a necessary condition for the initial value
problem y'=f(t,y) to have a solution on a given domain.
The second interesting question in handling first-order ODEs is to ask about the
uniqueness of solution. The guarantee for uniqueness of solutions becomes important,
especially when we want to find an approximate solution for the problem. Otherwise, our
computation may never converge. Shampine and Gordon [99] claim that if f(t,y) is
continuous and satisfies a Lipschitz condition on an open regiOn
D={(t,y) I a ::; t ::; b, -oo < y < oo }, then the problem y'=f(t,y), y(a)=A has a unique
solution for all intervals [a,b]. Braun [10] used the same premises as those in the
existence solution, in order to show the uniqueness of a solution of a first-order ODE.
1.3 Taylor's Series
Newton and Taylor has introduced a monumental concept of approximating a
function using a polynomial generated involving the derivatives of the function. This
5
method of approximation has been used in many applications. Unfortunately, in order to
gain a cheap cost of computation, this method is rarely used by itself in solving problems.
But, some methods in ODE computations such as Euler, Runge-Kutta, etc. use this
Taylor's series expansion. By some mathematical manipulations, the need to calculate
derivatives occurring in Taylor's series can be omitted. The expansion of a function f
using Taylor's n-series needs the condition that f should have n+ 1 continuous derivatives.
Suppose f(t) has n+1 continuous derivatives in [a,b]. For a point t, and c in [a,b], the
Taylor expansion off(t) around cis given as
f(t) = f(c) + f'(c)(t-c) + f 2>(c)(t-c)2/2! + ... + rn>(c)(t-ct/n! + Rn+I(~),
where Rn+I(~) = rn+I)(t-ct+1 I (n+1)!, and ~ is a point between t and c. Taylor's
expansion of function f of two variables t and y, where f and all its partial derivatives of
order up to n+ 1 are continuous in a neighborhood of the point ( a,b) is given as
f(t,y) = f(a,b)+ I ar+sf(a,b) ct-af (y- b/
~r s +R
l:5r+s:5n Ul 8y r! S! n+l'
where
R = " 8r+sf(~,TJ) (t-a)r (y- b/
n+l L..,. r s
r+s=n+l Ot 8y r! s! '
with the value (~,11) situated on the line segment joining the points (a,b) and (t,y).
Barton [8] has used this method to solve stiff problems. He equipped the method
with five devices: a device for estimating the local error, a device for predicting a stepsize
h, the use of predictor and corrector formulas, a device for deciding the order of Taylor
series to be used, and a device to detect when a transient has no effect on the solution.
6
Corliss and Chang [24] have developed software, ATSMCC (Automatic Taylor Series by
Morriss, Chang, and Corliss), using this Taylor approach. The software can be used to
solve nonstiff systems and moderately stiff systems of initial value problems of ODEs.
1.4 Stability and Instability of ODEs
All numerical methods are proposed to approximate the true solutions of
mathematical problems that are not easy to solve. In order to have accurate solutions, the
methods should be stable, so we can approximate the true values by controlling the errors.
When trying to solve an ODE problem, sometimes we are facing with instabilities of
solutions. The instabilities are not only caused by the approximation method of solution,
but are caused by the ODE problem itself. In general, we can classify instabilities of
ODEs solutions into two categories: inherent instability and induced instability.
Inherent instability is an instability caused merely by the property of the
differential system, and not the method of solution, while induced instability is an
instability caused by the use of a particular numerical method in solving the problem.
We will discuss the latter in the next chapter. To have a better understanding of inherent
instability, let's take a look at the following ODE problem: y" - y' - 6y = 0, with initial
conditions y(O) = 3, and y'(O) = -6. The analytical solution of this problem is y(t) = 3e-21
.
This solution will tend to zero when t goes to infinity. Suppose we change the initial
value y(O) by adding to it a small value f: > 0, so that y(O) = 3 + E. The analytical solution
for this new initial value is y(t) = (3+ 3E:/5)e-21 + (2E/5) e31
• It can be shown that no matter
how small f: is chosen, the solution will tend to infinity as t goes to infinity. In this case,
7
the solution y(t) = e-21 is said to be unstable. The picture of the solution forE = 0, 0.01,
and 0.1 are given in Figure 1.1. In numerical analysis, the situation when a little change
in the initial condition gives an unstable solution is called an ill-conditioned problem. It
is difficult to solve an ill-conditioned ODE problem numerically since the numerical
method usually has rounding errors. The approximate solution will always tend to
infinity just like the example of changing an initial condition given above. It is important
then to look at the stability of the system before performing some numerical
computations.
10
8
I
6 ....j
Vl ·x
<( I >-
4
2
0
0.0
y"-y'-6y=0, y(0)=3+E, y'(0)=-6
• E = 0 .
• E = 0.01
..... E = 0.1
----~
0.5 1.0 1.5
tAxis
Figure 1.1 Solution y = (3+ 3E/5)e-21 + 2E/5e31
forE= 0, 0.01, and 0.1
,.
• ___.•
2.0
8
The concept of stability of a system given here is based on the concept written by
Braun [1 0]. The question whether the solution ~(t) of the differential equation y'=f(y),
with y(O) = y
0
is stable or unstable is just the same as asking whether every solution \jl(t)
of y '= f(y) starting from a value that is very close to ~(t) at t=O remains close to ~(t) for
all t in the given interval or not. The solution y=~(t) of y ' = f(y) is said to be stable if
every solution \jl(t) of y ' = f(y) which starts sufficiently close to ~(t) at t=O must remain
close to ~(t) for all future times t. The solution ~(t) is unstable if there exists at least one
solution \jl(t) ofy' = f(y) which starts near ~(t) at t=O but which does not remain close to
~(t) for all future time. In his book, Braun said that if A is a matrix of ODEs y ' = Ay then
1. Every solution y = ~(t) of y ' = Ay is stable if all the eigenvalues of A have
negative real part.
2. Every solution y = ~(t) of y ' = Ay is unstable if at least one eigenvalue of A
has positive real part, and
3. Suppose that all the eigenvalues of A have real part ~ 0 and A-1 = icr 1, . • .,
Am= icrm have zero real part. Let A.j = icrj have multiplicity kj. This means that
the characteristic polynomial of A can be factored into the form
p(A.) = (A.-icr
1
)k 1
••• (A.-icrm)kmq(A.) where all the roots of q(A.) have
negative real parts. Then, every solution y = ~(t) of y ' = Ay is stable if A has
kj linearly independent eigenvectors for each eigenvalue A.j = icrj. Otherwise,
every solution ~(t) is unstable.
9
1.5 Convergence of Euler Method
Euler method is the simplest approximation method for the initial value problem
y' = f(t,y) with y(t0) = y0. This method is called the simplest method because it is easy to
derive directly from Taylor's series or directly from integration of the function by
assuming the tangent f(t,y) is constant at each interval. In practice, this method is not
recommended. One reason for not recommending this method is because this method has
a lack of accuracy compared to other methods at the equivalent stepsize. If f(t,y) is
assumed constant over a small interval (ti>ti+I) , a simple numerical procedure can be
derived for calculating numerical values YI> y2, .•••• ,yn, that approximate the true solution
values y(t1), y(t2), •.• , y(tn), ofy' = f(t,y) respectively. By integrating the constant value
ti+l
f(ti,yi) over the interval (ti>ti+1), we have y(ti+1) = y(t)+ Jf(ti,y(t))dt. The result of this
I;
integration gives the relation Yi+I = Yi + f(ti,yi)(ti+I-ti) where Yo= y(t0), which is Euler
method. If we examine this relation, it is none other than the first two terms of a Taylor's
series of y(t). A graphical interpretation of this method is shown in Figure 1.2, where
y ' =f(t,y) is interpreted as the slope of the integral curve at any point (t,y) on the curve.
Here, the solution curve on the range of (ti,ti+I) is approximated as a straight-line which
passes through (ti,yi) with the slope f(ti,yJ This type of explicit Euler method is
sometimes called "the crude Euler method".
The crude Euler method can be corrected by choosing the slope of the curve on
the range of (ti>ti+I) as the average of the slopes of the curve y(t) at the range endpoints.
The result of integration over the range of (ti>ti+I) then will give the approximate solution
10
of y ' = f(t,y) as Yi+I = Yi + h/2[f(ti,yi) + f(ti+I,yi+I)]. Since Yi+I appears implicitly on the
right-hand-side, this representation is called an implicit representation. This corrected
Euler method itself is sometimes called the Trapezoidal method. The value of Yi+I
appeared in f(ti+I,yi+I) is calculated using the crude Euler method. The procedure of
finding the value ofyi+I at ti+I by the Trapezoidal (corrected Euler) method is given as:
Predict Yi+I using the crude Euler method, y~!~ = Yi + hf(ti ,yJ.
h
Correct Yi+I using the Trapezoidal method, y~=~ = yi + 2[f(tpy) + f(ti+I,y~:~)].
y(ti) •················································································
Yi+I
Yi ···················································
slope=f(ti,yi)
ti ti+l t
Figure 1.2 Crude Euler Method
The formula for predicting Yi+I is usually called a predictor, while the implicit
.L...... •. _
11
formula for correcting Yi+I is called a corrector. The procedure of predicting a result by
one formula and then correcting by another formula is one of the most effective methods
for solving initial value problems and is called a predictor-corrector method. To get a
better result, the latter procedure can be used repeatedly until \Y~:~ - y~!~ \ <E. Let us
take y ' = 4.e0
·
81-0.5y with y(O) = 2 as an example problem solved both by the crude Euler
and the Trapezoidal methods. The exact solution for this problem is
y(t) = (4eo.st- 1.4e-o.st)/1.3.
Without need to iterate the Trapezoidal method, its result is still better than that of the
crude Euler method. Figure 1.3 is a graphical results of computations for stepsize
h = 0.5. It can be seen that the Trapezoidal method is better than the crude Euler method.
35
30
25
20
~"'
>-
1 5
1 0
5
0
0.0
y'=4.e 081-0.5, y(0)=2
• Exact 1
/,.
/
/
P~// • Euler ..
#
~/
... Trapezoidal
~
0.5 1 .0 1 .5 2.0 2.5 3.0
tAxis
Figure 1.3 Numerical Results Using
Exact, Euler, and Trapezoidal Methods
12
In order to have an acceptable computation, the analysis of errors caused by the
numerical solution of ODEs becomes important [66, 78, 80]. Usually the errors are
classified into two types of errors: round-off errors and truncation or discretization errors.
Round-off errors are errors that occurred because of the limited numbers of significant
digits that can be retained by a computer. Truncation errors consist of a local truncation
error that results from computations over a single step, and a propagated truncation error
that results from the previous steps. The sum of the local truncation error and the
propagated truncation error is called a total or global truncation error [15, 44, 57, 66, 80].
Golub and Ortega [49] have shown that if the function fhas a bounded partial derivative
with respect to its second variable, and if a solution y 1 = f(t,y) with y(t0) = Yo has a
bounded second derivative, then the Euler approximation converges to the exact solution
as h ~ 0, and the global discretization error of Euler method satisfies E(h) = O(h).
1.6 The Multistep Concept
The numerical methods of solving an ODE y 1 = f(t,y) with initial value y(t0)=y0
are classified into two categories: single step methods and multistep methods. In single
step methods, the approximate value ofy(t) at ti+l is calculated merely using the values Yi,
y/ and h. In multistep methods, the approximate value of y(t) at ti+l is calculated using
the recurrence relation in terms of the function values y(t) and derivatives Y1(t) at ti+l and
at previous nodal points. Among examples of single step methods [66, 77, 80] are:
1. Yi+l = Yi + f(ti,Yi) (Euler method)
13
y: Y(n)
2. Y· I= Y· +-1 h+ + _I_hn I+ I 1.I ... n!
(Taylor's series method)
3. Yi+t = Yi + hf(ti+h/2,yi+(h/2)f(tbYi)) (Midpoint rule)
Among examples of multistep methods are:
h
1. Yi+l =yi +2[3f(ti,yJ-f(ti-PYi-1)] (Adams-Bashforth)
h
2. Yi+l =yi +2[f(ti,yJ+f(ti+PYi+l)] (Adams-Moulton)
The idea of multistep methods appears when some information from previous
points has been gained using single step methods. After gaining the information, the
value of f(t,y) at the next step is based on interpolation over that information.
Computation using these multistep methods is in general more accurate than those in
single step methods. Suppose we already have information of y1, y2, ... , Yn using single
step methods. By that information, we can calculate f0, f~> .... , f"" A multistep method
tn+l
of solving y' = f(t,y) can be based on the relation y(tn+l) = y(tn_M) + Jf(t,y)dt, where M
tn-M
is some non-negative integer. In this case we will approximate the function f(t,y) by a
polynomial P(t) generated using information f"' fn-~> ... , fn-M· Later on we will discuss
another way of finding multistep methods. That is, instead of approximating the
integrand f(t,y), the methods will be achieved by approximating y(t) using a polynomial
and then by differentiating the polynomial. By substituting this interpolation polynomial
tn+l
P(t) for f(t,y) as the integrand, we have the relation y(tn+l) = y(tn_M) + JP(t)dt. In order
tn-M
14
to integrate this from tn-M up to 1n+I we should extrapolate P(t) through 1n+I· The
procedure of solving y' = f(t,y) is then similar with that in the Trapezoidal method:
tn+l
1. Predict Yn+l using y:
1
: 1 = Yn-M + JP(t)dt.
1n-M
2. Compute fn+l using fn+l = f(tn+P Y~~~) ·
3. Interpolate P(t) using fn+I> f"' fn-1> ..... .
tn+l
4. Correct Yn+l using y::: = y n-M + JP(t)dt.
tn-M
This procedure can be executed repeatedly until we satisfy with the desired accuracy of
the corrected result. Suppose E > 0 is the accuracy needed for the computation. Then
steps 1 through 4 will be executed repeatedly until IY~:~ - y~~~ I< E. Usually, P(t) is
approximated using Newton's backward-form (NBF). The followings are some examples
of predictors and the related correctors [66, 78, 80]:
Milne method:
Predictor old 4h
: Y n+l = Y n-3 + 3 [2fn - fn-1 + 2fn-2] ·
Corrector
h
: Y~:~ = Yn-1 +3[fn+l +4fn +fn-1].
Adams-Moulton methods:
1. Second-order:
Predictor old h
: Yn+l =Yn +2[3fn -fn-1].
Corrector new h
: Y n+l = Y n + 2[fn+l + fn].
15
2. Third-order:
Predictor
old h
: Yn+l = Yn +U[23fn -16fn-l +5fn_2 ].
Corrector new h
: Yn+l =Yn +U[5fn+l +8fn -fn-1].
3. Fourth-order:
Predictor
: Yon l+dl = Yn + h
24
[55fn- 59fn-l + 37fn_2 - 9fn_3 ].
Corrector
h
: Y::7 = Y n + 24 [9fn+l + 19fn - 5fn-l + fn-2].
4. Fifth-order:
Predictor
old h
: Yn+l = Yn +
720
[190lfn -2984fn-l +2616fn_2
-1274fn-3 + 251fn-4 .
Corrector
new h
: Y n+l = Y n +
720
[25lfn+l + 646fn - 264fn-l
+ 106fn-2 -19fn-3l·
5. Sixth-order:
Predictor
old h
: Y n+l = Y n +
1440
[ 4277fn - 7923fn-l + 9982fn_2
-7298fn-3 + 2877fn-4- 475fn-sl
Corrector
new h
: Yn+l = Yn +
1440
[475fn+l + 1427fn -798fn-l
+ 482fn-2 -173fn-3 + 27fn-4l·
__.._____
CHAPTER II
STABILITY ANALYSIS
In the previous chapter the instability caused by the problems themselves has
been addressed. The next observations will be about the instability caused by the
approximation methods used to solve the problem. This instability is called induced
instability. As mentioned in the previous chapter, before using the approximation
method, we have to make sure that the problem does not belong to the class of illconditioned
problems. Due to the inaccuracy of the numerical methods used in solving
an initial value problem for an ordinary differential equation, errors occur at each step of
the integration. The total error then becomes the summation of the local truncation errors
and their propagation. If the method is unstable, the accumulation of these local errors
will make the total error becomes larger than it would be expected. This is the reason
why before using the approximation methods, we need to make sure that the problem
does not have inherent instability. The phenomenon of the growth of this error is called
numerically induced instability.
The stability of numerical methods of solving ODEs can be analyzed using the
simple linear first order differential equation proposed by Dahlquist [27, 28, 29]
Y' = 'Ay, y(to) =Yo (2.1)
16
17
where A is a constant. This ODE is also called a test problem which was first introduced
by Dahlquist [27, 28, 29] and is sometimes called Dahlquist's test problem. The solution
of (2.1) is y(t) = y 0e).(t-to). If tn is written as tn = t0 + nh, y(t) can be given as
y(tn+I) = e"hy(tn). The value of e"h can be approximated using Taylor's series. If the
approximate value lh is denoted by E(Ah), then y(tn+I) can be approximated by
Yn+I = E(Ah)Yn· Assume E(Ah) is taken as a Taylor's series of order p, then
E(Ah) = 1 + Ah + (Ah)
2
+ (Ah)P
2
' ... + ---
. p!
(2.2)
Based on the approximate solution given above, it can be said that a numerical
method is stable, if the error Yn-y(tn) = en remains bounded as n~oo. By calculating the
error resulted from the approximate solution, we have
en+l =yn+l -y(tn+l)
= E(Ah)[y(t") +en]- e'.hy(t") (2.3)
= [E(Ah)- e'·h ]y(tn) + E(Ah)en.
From (2.3), it can be shown that the error at tn+I consists of two parts. The first part,
E(Ah)- e'·h, is the local truncation error that can be chosen as small as we want by
choosing a suitable approximation of E(Ah). The second part, E(Ah)en is the propagation
error from the previous step tn to tn+I that will not grow if \E(Ah)\ ~ 1.
Definition 2.1
For a given hand A from Dahlquist's test problemy ' = Ay, a numerical method
Yn+I= E(Ah)Yn is called absolutely stable if\E(Ah)\ ~ 1.
18
Definition 2.2
For a given hand A from Dahlquist's test problemy 1 = Ay, a numerical method
y n+ 1 = E(Ah)y n is called relatively stable if IE(Ah)l ~ el.h.
The region of stability of numerical methods will be based on the methods used to
approximate the solutions [2, 4, 6, 14, 25, 66, 79].
In the case of an ODE in the form of y 1 = f(t,y), with initial value y(t0) = y0, the
problem can be transformed using the Taylor expansion about (tn,Yn), so that the problem
can be solved as a linear problem. If the Taylor series is truncated after first order terms,
the linearized form of the equation will be given as:
y I= Ay + Bt + c (2.4)
where
or or af af
A=8-y( t n' yn) ' B=8-t (t n' yn) ' C=fn -y n -8y( t n' y n )-t nO-t( t n' yn) ·
The solution of (2.4) can be found in the form y = Ae'·t + p(t), where p(t) is a
polynomial of degree at most 2. It can be seen that the solution is dominated by the
exponential term.
2.1 Stability of Euler Methods
The crude Euler method is Yn+l = Yn + hf(tn,Yn). If we apply the test problem
y 1 = Ay, we will have Yn+l = (l+Ah)Yn- In this case E(Ah) = 1+ Ah. It can be said then that
the crude Euler method is absolutely stable if 11 +Ahl ~ 1. For A real, the interval of
___j,____ -~
19
absolute stability is -2 ~ Ah ~ 0. For A a complex number, the region of stability is the
region inside a circle with center ( -1 ,0) and radius 1 as shown in Figure 2.1.
The implicit Euler method is Yn+I = Yn + hf(!n+t>Yn+I)· If we apply Dahlquist's test
problem y ' = Ay, we will have Yn+l = Yn + hAYn+l or Yn+l = (1-Ah)"1Yn· In this case
E(Ah) = (1-Ah)"1
• It can be said then that the implicit Euler method is absolutely stable if
11-Ahl ~ 1. For A real, the interval of absolute stability is {Ah I Ah ~ 2 or Ah ~ 0}. For A
complex, the region of stability is the region outside the circle with center (1,0) and
radius 1 as shown in Figure 2.2. To have a better understanding of this stability concept,
let takes an example y ' = -50y, y(O) = 1. For the implicit Euler method, the numerical
computation can be done by applying a predictor function and then applying Newton's
iteration.
Im(A.h)
II+hA.I ~l
-2 Re(A.h)
Figure 2.1 Stability Region of the
Explicit Euler Method
Im(A.h)
ll-A.h I~ 1
2 Re(A.h)
Figure 2.2 Stability Region of the
Implicit Euler Method
20
The algorithm of the implicit Euler method using a predictor-corrector scheme is as
follows [46, 72, 94]:
1. Predict Yn+I using the crude Euler method, y~~1 = y" + hf(t", y").
2 . C orrectyn+I us.m gth e1.m p11"C .1 tmeth od, Ynn+ewl =yn + hfi1~ _tn+l'Yno+ldl ) ·
The algorithm ofthe implicit Euler method applying Newton's iteration is as follows [51,
70, 72]:
1. Give the initial value for y n+ b says y gues, so that Y~1~1 = Y gues ·
old old old
2. Compute F(Yn+J), F(y n+l) = Y n+l - Y n- hf(tn+l 'Y n+l ).
3. ComputeF'(Yn+I), F'(y~1~1 )=1-hf, (tn+l'y~1~1 ) •
., n+l
4 C
· N , · · new old F(y~~l)
. ompute Yn+l usmg ewton s IteratiOn, y m+l = y n+l - T"'l' old
Y n+l
5. Iterate again from step 2 untiljy~:'; - y~~1 j < EPS, where EPS is an accuracy
requested to a numerical solution.
The exact solution of this problem is y = e-sot. A graphical representation of the explicit
Euler method for h = 0.01, 0.03, 0.04, and 0.05 is given in Figure 2.3. The computation
of y ' = -SOy, using the Euler method gives results as follows :
- The solution converges to the exact solution for every stepsize which is less
than 0.04. If the stepsize is close to 0.04, the solution will oscillate but it will
converge to the exact solution.
- For stepsize h=0.04, the solution will never converge to the exact solution
even though the computation is stable. The solution just oscillates between
21
minus and plus the initial value by turns.
For every stepsize larger than 0.04, the computation will never converge and
it will oscillate and then go to minus or plus infinity by turns.
8'
6
4 r 0 Exact
v h=O.Ol
c:: h=0.03
~ h=0.04
& h=O.OS
2
y'=-50y, y(O)=l
·(~x/ o)o.l .~bo~ .--Q , ~ ~~ o ~~ o o o ¥o l)-v~
!
-2'
I
i
I
-4~
-6
i -8 L
Figure 2.3 Euler Method with Stepsize
h=O.Ol, 0.03, 0.04, and 0.05
22
The comparison of computations among crude Euler, implicit Euler usmg
predictor-corrector, and implicit Euler using Newton's iteration has been made for the
case of problem y 1 = -SOy, y(O) = 1. The conclusions taken from the computational
experiments for the case y 1 = -SOy are:
1. The implicit Euler method using Newton's iteration always converges to the
exact solution.
2. As mentioned in the previous experiments, the crude Euler method always
converges to the exact solution for a timestep less than 2/SO = 0.04. For the
timestep h=0.04, the method is still stable, but does not converge to the right
solution.
3. The computation of the implicit Euler method as a corrector and crude Euler
as predictor is convergence to the exact solution for a timestep less than
1/SO = 0.02 (see Appendix E). For the timestep h = 0.02, the computation is
still stable, but the result is not convergence to the right solution. In the crude
Euler method the computation still converges to the exact solution for
timestep h = 0.03 (even though oscillating), but the computation using this
kind of predictor-corrector does not converge for timestep h ~ 0.03. It just
goes to infinity.
If the computation using the implicit Euler method as corrector and the crude Euler as
predictor is stable, then the result converges to the result given from the implicit Euler
using Newton's iteration. This PC method may not converge yet for one iteration, but if
23
we repeatedly correct the computation, the result will be the same with the result from
the implicit Euler using Newton's iteration.
Graphical representation of the comparison among the three methods aforementioned for
stepsizes h = 0.01, 0.03, 0.04, 0.05 are given in Figure 2.4.
J
1.00 20
18
/
h=O.OJ
c Exact
" Euler I 0 Euler+PC
t:. Euler+Newton
H /
:~ ~ 0.~0 0.05 0.10
-2 L
16
0.75 14 ll\ h=O.OI
Ul
.~ l\\\ ')(
11) i :>. 0.50 i 0 E:uct
? Euler
c Euler+ PC
6 Euler+Newton
I 1\\\
0.25
12
10
8
6
~
0.15 0.20 0.25
0.00 I ~·~
t-axis
0.00 0.05 0.10 0.15 0.20 0.25
t-axis
20
16 ~
I
16
12 ~
h-0.04 12 ~ I b=O.OS
I I c Es.act
Exact
v Euler
Ul I " Euler 8 ~ 0 Euler+ PC
')(
8 "
c Euler+ PC
Ul I 6 Euler+Newton
11)
I
6 Euler+Newton
')(
I >-
11)
I i >- 4 ~
4 L
I I Up'~ u ~.;y 0.;0 \ 0,;
t-ax~ -4 ~
.30
t-axis -8 ~
Figure 2.4 Crude Euler, Implicit Euler Using PC,
and Implicit Euler Using Newton's
24
Graphical computations ofy'=-50y, y(O)=l using implicit Euler with Newton's iteration
for h =0.01, 0.04, and 0.07 are given in Figure 2.5, where the gradient is also shown.
y-ax1s
1
0.8
0.6
0.4
0.2
y•=-50y, y(0)=1
~ ~ f f f f f f f ~ f f f ~ f f f t f f
\ \ \ ,\ ~~ \ \ \ ~ \ \ \ \ \ \ \ ~ \
~,\ ~ \ ff\J.Ut \ \ \ ~ \ \ \ ~ ~ \ \ ~ ~
\ \ \ \ \ \ \ ~ \ \ \ ~ \ \ \ ~ ~
\ \ \ ~ \ \ \ \ \ \ \ \ \ ~ \ \ \ ~ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
0.05
"t,.....)) {\1 \ \, \, \, \ \, \, \, \, \, \,
-trv·'-{, .... ... ... ... ' ... ... ... \ ... ...
\ ' ~ ~ \ ' ~ ~ \ \ ~ ~
~ h=O.Q4 ~ ~ \ \ ~ ~
~ ~ ' ' ' ~ ~
... ... ... ... ... .. .. .. .. ..
0.1 0.15 0.2
Figure 2.5 Implicit Euler by Showing the Vector Gradient
2.2 Stability of Explicit Runge-Kutta Methods
t-axis
When Carl Runge and Wilhelm Kutta introduced Runge-Kutta methods, they
were motivated by Euler method and Taylor's series. There are three major
characteristics of Runge-Kutta methods:
25
- Single-step methods,
- Derived from Taylor's series with the same order, and
- A voiding the evaluation of derivatives of f.
While the Euler method only evaluates the function f(t,y) of y ' = f(t,y) at the solution
points themselves, the Runge-Kutta method is done by evaluating the function f(t,y) at
several points within the integration interval. Let us consider the second-order Runge-
Kutta method that is given as Yn+I = Yn + ak1 + bk2, where k1 = hf(tn,Yn),
k2 = hf(tn +Ah,yn +Pk1). The constants a, b, A, and P are determined so that the equation
will agree with the Taylor expansion. The Taylor expansion ofy(t) at 1n+I through terms
of order h3 is
fi 2 3
Y(1n+I) = Y(1n) + h (1n,Yn) + h /2[ft(1n, Yn) + f(1n,Yn)fy(1n,Yn)] + O(h ).
By doing the same thing to f(t,y) we get
f(tn +Ah,Yn +Pk1) = f(1n,Yn) + Ah~ (1n,Yn) + Pklfy (1n,Yn) + (A
2
h
2
/2)fu (1n,Yn) +
2 2 3 APhkifty(1n,Yn) + (p k1 /2)fyy (1n,Yn) + O(h ).
By substituting the latter expansion into Yn+1=yn+ ak1 + bk2, we have
2 3
Yn+l = Yn + (a+b)hf(1n,Yn) + h [bAftC1n,Yn) + bpf(1n,Yn)fy (1n,Yn)] + O(h ).
This latter equation will be identical with the Taylor expansion of y(t) at tn+I through
terms in h2 if a+b = 1, bA = 1/2, and bp = 112. If a =112 then b = 1/2, A= 1, p = 1, and
Y n+ 1 is written as
Yn+l = Yn + 112[kl + k2] (2.5)
where k1 = hf(tn,yn) and k2 = hf(tn +h,yn +k1). Equation (2.5) is called a second-order
Runge-Kutta method. Its graphical interpretation is shown in Figure 2.6. The stability of
26
second-order Runge-Kutta method can be examined using Dahlquist's test problem
y' = 'Ay. By substituting the test function into (2.5), we have
2
Yn+l = Yn + l/2[h'Ayn + h'A(Yn + h'Ayn)] = [1 + 'Ah + ('Ah) /2]Yn ·
Based on Definition 2.1, it can be said that the second-order Runge-Kutta method will be
stable if Jl + 'Ah + ('Ah)212J ::; 1. A graphical representation of the region of stability when
h'A is complex is given in Figure 2.7. For h'A a real number, the interval of stability is
-2 ::; h'A ::;o.
k2
/ -?
Yn+l .]......................................... __;;,
Yn 112[kl+k2]
kl
ln ln+I
Figure 2.6 A Second Order ofRunge-Kutta
If a is chosen 1/4, the method is called the Heun method and is written as
Yn+l = Yn + kl/4 + 3k2/4 (2.6)
where k1 = hf(tn,Yn), and k2 = hf(tn+2h/3,yn+2k1/3).
With the same approach, Lapidus and Seinfeld [80] have showr1 :
1. Third-order Runge-Kutta
Yn+l = Yn + 1/6 [kl + 4k2 + k3] (2.7)
.....
27
where k1 = hf(tn,Y0
), k2 = hf(tn+hl2,y0+k1/2), and k3 = hf(tn+h,y0-k1+2k2). The region of
absolute stability for A.h complex is
{ A.h Ill+ A.h + (A.h)
2
12 +(A.h)
3
161 ~ 1}
(see Figure B-1 of Appendix B). For A.h a real number, the interval of absolute stability
satisfies -2.512745 ~ A.h ~ 0.
Order-2 Im(A.h)
Figure 2.7 Stability Region of Second Order
Runge-Kutta
2. Fourth-order Runge-Kutta
Yn+l = Yn + 1/6 [kl + 2k2 + 2k3 + k4] (2.8)
where k1 = hf(tn,Y0), k2 = hf(tn+hl2,yn+kl/2), k3 = hf(tn+hl2,yn+k2/2), k4 = hf(tn+h,yn+k3).
The region of absolute stability for A.h complex is
{A.h Ill+ A.h + (A.hil2 +(A.h)3/6 + (A.h)
4
1241 ~ 1}
(see Figure B-1 of Appendix B). For A.h a real number, the interval of absolute stability
satisfies -2.78529 ~ A.h ~ 0. For A.h real numbers, the intervals of absolute stabilities of
......
28
explicit Runge-Kutta methods for orders 1 up to 10 are shown in Table 2.1. Butcher [15]
proposed an algorithm to generate data for plotting the stability region of Runge-Kutta
methods. We developed FORTRAN (see Appendix B) and Mathematica (see Appendix
F under the name of "Mathematica 5") codes with a different approach. The approach is
based on the implementation of Newton-Raphson method for finding a root of a
polynomial. Here, the polynomial is evaluated using nested multiplication which is
recognized as Homer's rule.
TABLE2.1
Stability Region ofRunge-Kutta
for Orders 1 to 10.
Order Interval for A.h real numbers
1 SL -2 ~ A.h ~ 0
2nu -2 ~ A.h ~ 0
3m -2.5127453266 ~ A.h ~ 0
4m -2.78529 ~ A.h ~ 0
5m -3.21705 ~ A.h ~ 0
6Ln -3.55344 ~ A.h ~ 0
t" -3.95413 ~ A.h ~ 0
8m . -4.31363 ~ A.h ~ 0
9Ln -4.70083 ~ A.h ~ 0
lOrn -5.06952 ~ A.h ~ 0
I
......
29
The general expression of as-stage Runge-Kutta method for solving y' = f(t,y),
y(t0) = t0 is given as:
s
Yn+l = Yn + Lwiki
i=l
i-1
wherek1 =hf(tmYn)and ki =hf(t 0 +cih,yn + Laijk), i=2,3, ... , s
j=l
(2.9)
To maintain consistency, which will be discussed in the next chapter, these
following conditions should be satisfied:
i-1 s
ci = Laij and Lwk = 1 , i = 2,3, ... ,s (2.1 0)
j=l k=l
Here, the w' s represent the weighting coefficients on the slopes, s is the total number of
stages, the c' s represent the subinterval locations at which the derivatives are being
evaluated, and the a's represent the slope weightings for the intermediate stages.
Butcher proposed a condensed representation of the Runge-Kutta method as
shown in Table 2.2.
TABLE2.2
Butcher's Array ofRunge-Kutta Methods
0
ci I a21
c2 a31 a32
Cs ~I as2 ass-!
WI w2 Ws
30
This representation is recognized as Butcher's array. The Runge-Kutta (2.8) method is
represented in Butcher's array as shown in Table 2.3.
TABLE2.3
Butcher's Array of a Fourth-Order Runge-Kutta Method
0
1/2 I 1/2
1/2 0 1/2
1 0 0 1
116 113 113 116
2.3 Stability of Implicit Runge-Kutta Methods
The general expression of a s-stage implicit Runge-Kutta method for solving
y' = f(t,y), y(t0) = t0 is given in Equation (2.11) as follows:
where
and
s
Yn+l = Yn + Lwiki
i=l
s
ki = hf(t 0 + cih, y n + Laijk j), i = 1,2,3, ... , s
j=l
s s
ci = Laij and L wk = 1 , i = 1,2,3, ... ,s.
j=l k=l
(2.11)
31
Some derivations of implicit Runge-Kutta methods have been shown by Jain [66]. A
second-order implicit Runge-Kutta method is given as
Yn+l =yn+kl (2.12)
where k1 = hf(tn+h/2,yn+k1/2). If we apply the test problem y ' = Ay, we have
k1 = h"J....(yn+k1/2). After finding kl> its value is substituted into Eq. (2.12) to get
l+"J....h/2
Y n+ 1 = 1 - Ah / 2 Y n . (2.13)
For Ah real, the interval of absolute stability satisfies J....h ::; 0, and for Ah complex, the
region of absolute stability satisfies I (1 + J....h/2)/( 1- J....h/2)1 ::; 1. Using a similar approach,
a third-order implicit Runge-Kutta method can be written as
Yn+l = Yn + 1/4 [ 3kl + k2] (2.14)
with k1 = hf(tn +h/3,yn +k1/3) and k2 = hf(tn +h,yn +k1). For Ah real, the interval of absolute
stability satisfies J....h ::; 0, and for J....h complex, the region of absolute stability satisfies
1[1 +2 J....h/3+( J....h)2/6]/[1- J....h/3]1::; 1.
2.4 Difference Equations
If y is a function oft and t is replaced by tk = t0 + kh, where k is an element of
{ .. ,-2, -1,0,1,2, ... }, then the dependent variable y(t) can be replaced by Yk as an
approximation for y(tk). The forward difference of the Yk value is denoted by L1yk and
defined as Yk+I-Yk· This definition implies that
2 L1 Yk = L1(L1yk) = L1(Yk+I-Yk) = Yk+2- 2yk+I + Yk= L1Yk+I-L1Yk·
The first derivatives of y overt are defined as
y(t +h)- y(t)
h
dy =lim
dt h--+0
y(t)- y(t- h)
h
y(t +h)- y(t- h)
2h
The finite difference approximation of the first derivative of y is defined using the first
two terms ofthe Taylor series ofy as
dy ~
dt
Yk+l -yk
h
Yk- Yk-1
h
Yk+l -yk-1
2h
forward Euler.
backward Euler
central difference scheme_
The interpretation of these methods in solving a problemy'= f(y) is as follows:
y k+l - y k = f(y k)
h
Yk- Yk-l = f(yk)
h
y k+l - y k-l = f(y k)
2h
forward Euler
backward Euler
central difference scheme
Let's take an example of the decay equation dy/dt = -y, y(O) = y0. The correspondence
of methods is
Yk+l- Yk = -yk
h
forward Euler
Yk -yk-l =-yk
h
backward Euler
= -yk central difference scheme
33
The exact solution of this problem is y(t) = y0e-t. The absolute value of the exact
solution monotonically decreases to zero as t ~ oo. The forward Euler scheme for this
problem is Yk+l = (1-h)Yk· By substituting the value ofy0, we have the relation Yk = (1-
h)k y0. The backward scheme for this problem is Yk = Yk-l - hyk or Yk = [1/(1 +h)] Yk-1·
By taking k = k+1, we can have Yk+l = [1/(1+h)]Yk· By substituting the value ofy0,
we have the relation Yk = [1/(1 +h)t y0. The relation for the central difference scheme is
Yk+l = Yk-l- 2hyk. The value ofy1 for this method is calculated using the Euler method.
All computations of this problem were done using the three methods, each for
h = 0.05,1.,2., and 3. From the observations shown in Figure 2.8, the backward method
is always stable for every 0 < h < oo. The central method is always unstable for every
h>O. For h = 1 and k>O, the value Yk of the forward method is identically equal to zero.
If 1 <h<2, the value Yk of the forward method decreases to zero with oscillating (change
in sign) amplitude. For h ~ 2, the Yk of the forward method ofthis problem does not lead
the solution to zero when t goes to infinity. The forward methods is always stable for h
in the interval 0 < h < 2. By applying the forward method to the test problemy ' = A.y
we have Yk+l = Yk +hA.yk = (1 +A.h)Yk· From the latter equation, for A.h real, the region of
stability of the forward method is an interval {A.h Ill +A.hl s 1} = -2 s A.h s 0 (similar to
the crude Euler method). In a similar way, the region of absolute stability of the
backward method for A.h real is an interval { A.h I A.h ~ 0 or A.h s 0} (similar to the
implicit Euler method). For the central method, the formula for finding the region of
stability is not simple.
1.
4
11=0.05 n ~ - ....,___..,;
0 Exact
\1 Forwanl (/) ·x (/) r \ r Backwanl
ctl
-~ 0.50
6 Central
>- . >-
\1 Forwanl
-sc- D Backwanl
I 6 Central
I
0.25' "Lt. -10 r
i
-12
-14
ooo I ~
0 1 2 3 4
taxis
15
12 b=J.
'~·:\·~ 9 c Exact
- \1 Forwanl
2
1- taxis
-4~ c Backwanl
i 6 Central
-8 0 Exact (/) 6 (/) \ ·x
\1 Forwanl ctl ·x
ctl -8 c Backwanl >- >- 6 Central
-10
·12
-14
-16
-18
-20
Figure 2.8 Difference Equation Methods Using Forward,
Backward and Central for h=O.OS, 1, 2, and 3
34
CHAPTER III
STABILITY OF MULTISTEP METHODS
We have seen that some ideas for solving Ordinary Differential Equations lead to
multistep methods. Unfortunately, the computations of multistep methods are not as
easy as those in single step methods, because in multistep methods some starting values
are needed but they are not given in the starting information. Traditionally, to support
this lack of information we can use Runge-Kutta methods with the same order accuracy.
There are two possibilities that can be used in order to provide starting values. First, use
a single step method to compute all starting values that are not given in the information.
Second, use a single step method to support information for a two-step method, then use
this two-step method to support information for a three-step method and so on until all
information needed to support the chosen multistep method is fulfilled. In general, the
multistep methods not only include explicit single step methods, but also implicit single
step methods. A general multistep or k-step method for the solution y ' =f(t,y), y(t0) =Yo
can be written as
aoy n + aly n+l + a2y n+2 + + ak.1Y n+k·I + aky n+k =
(3.1)
h~(tn 'tn+P .. 'tn+k·P tn+k 'Y~ 'Y~+ p ... 'Y~+k-P Y~+k; h)
where his the constant stepsize and ao, a~> .. , ak are real constants. If~ independent of
y n+k' then the general multistep method is called an explicit, open or predictor method;
35
36
otherwise it is an implicit, closed or corrector method. A general linear multistep or k-step
method can be written as
aoy n + aly n+l + a2Y n+2 + · · · + ak-1Y n+k-1 + aky n+k =
h(boY~ + b~y~+l + · · · + bk-IY~+k-1 + bky~+k)
(3.2)
This formula can be simplified to
p(E)y n - hcr(E)y n' = 0 (3.3)
2 k-1 k 2 k-1 k
where p(~) = ao + a1 ~ + a2~ + .. ak-I~ + ak~ , cr(~) = b0 + b 1 ~ + b2~ + .. bk-I~ + bk~ ,
and E is the shift operator which is defined as Ek Yn = Yn+k· Dahlquist [1956] was the
person who started to study a general theory of multistep methods and was the first
person who introduced the polynomials p(~) and cr(~).
In the same manner as single step methods, the theory of multistep methods also
involves the concepts of consistency, stability, and convergence. Dahlquist [27 ,28] has
defined that Equation (3.3) is said to be consistent if p(l)=O, and p ' (1) = cr(l). We
know that the initial value problemy ' = 0, y(O) = 1 has the exact solution y(t) = 1. If we
solve this problem using (3 .2), we have
aol + a11 + ... + akl = h(b00 + b10 + ... +bkO) or p(l) = 0.
Let's take another example, y ' = 1 with y(O)=O. The exact solution of this example is
y(t) = t. Substituting this exact solution into (3.2), we have
ao(nh) + a2(n+l)h + ... + ak(n+k)h = nh[a0 + a1+ ... + ak] + h[la1+2a2+ ...
+ (k-l)ak_1+kak]
= h(b0 + b1 + ... + bk)
...........____
37
where y 5 = sh, for s=O, 1, 2, .... n+k. Simplifying the latter equation and applying the
result of p(l)=O, we have p'(1) = cr(l). Dahlquist [27,28], Lapidus and Seinfeld [80],
and Jain [66] have shown that the multistep method (3.1) is of order p, if and only if one
of the following equivalent conditions is satisfied :
k k k
1. Lai = 0, and Laiiq = qLbiiq-I for q = 2, ... ,p,
i=O i=O i=O
2. p(eh)- hcr(eh) = O(hp+I) for h---+ 0, and
3. p(~) -cr(~) = 0((~-1)P) for~---+ 1.
log~
Definition 3.1
(3.4)
(3.5)
(3.6)
The formula (3.1) will be said to be absolutely stable in the sense of Dahlquist if
all roots ~i of the characteristic polynomial p(~)=O are such that l~d ~ 1 and those
roots for which l~d = 1 are simple.
Theorem 3.1
Consistency and stability are together necessary and sufficient conditions for the
formula (3 .1) to be convergence.
Proof: See Dahlquist [27,28].
Example 3.1:
Investigate the consistency and the order of a numerical method
Yn+2 = -4Yn+l + 5yn + h(4fn+l + 2f0 ).
Solution:
.....
38
For this case k = 2, ao = -5, a1 = 4, a2 = 1, b0 = 2, b1 = 4, b2 = 0, cr(s) = 2 + 4s, and
p(s) = -5 + 4s + s 2 • Since p(l) = -5 + (4)(1) + 12 = 0, and p'(l)= 4 + (2)(1) = cr(1), then
the method is consistent. From condition 2 (3.2), we have
p( eh) - hcr( eh) = -5 + 4eh + e2h - h(2 + 4eh)
= - 5 + 4[1 + h + h2/2 + O(h3
)] + [1 + (2h) + (2h)2/2 + 0((2h)
3
)
- 2h- 4h(l + h + h2/2 + O(h3
)]
= -5 + 4 + 4h + 2h2 + O(h3
) + 1 + 2h + 2h2 + O(h3
) - 2h - 4h - 4h2
-2h3
- O(h3
) = - 2h3 + O(h3
) = O(h\
So the numerical method is of order 2.
Example 3.2:
Investigate the consistency and the order of the trapezoidal method
Yn+l = Yn + (h/2)[fn + fn+l].
Solution:
In this case k = 1, ao = -1, a1 = 1, b0 = 112, b1 = 112, cr(s) = 1/2 + (1/2)s, and
p(s) = -1 + S· Since p(1) = -1 + (1)(1) = 0, and p'(l)=1= cr(1)=1/2+(1/2)(1), then the
method is consistent. Let us now use condition 1 to show the order of the trapezoidal
method. By applying condition 1(3.1), we can see that the condition will be satisfied
only for p = 2. So the numerical method is of order 2.
3.1 Linear Difference Equations
Among other purposes of studying difference equations are the use of these
equations to formulate and analyze discrete-time systems, to approximate the integrations
....
39
of differential equations using finite-difference schemes, and to study deterministic chaos
[91, 92]. In general, an ordinary difference equations is a relation ofthe form
Yn+k = F(n,yn,Yn+b · · ., Yn+k-1), (3.7)
where the order is the difference between the highest and lowest indices that appear in
the equation. Based on this definition, the order of Equation (3. 7) is k, and shifting the
labeling of the indices will not change the order of the equation. If Equation (3. 7) is
shifted to Yn+k+r = F(n+r,Yn+r,Yn+r+b ... ,yn+r+k-1), its order is still k. A solution of difference
equation (3.7) is a function <p(n) that reduces the equation to an identity. If we substitute
<p(n)= 2" into Yn+l - 2yn = 0, we will have an identity in the form of 2n+l - (2) 2" = 0. In
this case, <p(n)=2" is said a solution of Yn+l - 2yn = 0. A difference equation can be
derived from a sequence {Yn} to which is defined as a function of n and k arbitrary
constants Ct. c2, ... , ck. The difference equation of this sequence will be of order k. If
we calculate Yn+l from Yn = A2", we will have Yn+l = A2n+l = (2)A(2)" = 2Yn· The
difference equation ofyn = A2" is a first-order equation in the form ofyn+l - 2yn = 0. By
calculating Yn+l and Yn+2 from Yn = c12" + c22", we will have a second-order of difference
equation in the form of Yn+2 - 7Yn+l + lOyn = 0. This gives a difference equation
Yn =An+ f(A) is Yn = (Yn+l - Yn)n + f(Yn+l - Yn) which can be found by calculating Yn+l
and substituting the value of A back into the first equation.
A kth-order linear difference equation with constant coefficients is a relation in
the form of
Yn+k + a1Yn+k-l + a2Yn+k-2 + · · · + akYn = Rn (3.8)
40
where a!> a2, ... , and ak are constants with ak '* 0, and R0 is a function of n. The equation
(3.8) is called a kth-order nonhomogeneous linear difference equation with constant
coefficients. If R0 = 0, then Equation (3.8) is called homogeneous and it becomes
Yn+k + a1Yn+k-I + a2Yn+k-2 + · · · + akYn = 0. (3.9)
Using the shift operator E, we can write Equation (3.9) in the form f(E)y0 = 0, where
.c.(E) Ek Ek-I Ek-2 . h fi . Th h . . . 11 = + a1 + a2 + ... + ak 1st e operator unctwn. e c aractenst1c equatiOn
associated with Equation (3.9) is f(r) = 0.
Theorem 3.2
L et ri b e any so1 u t.w n toth e c h aracten.s t1.c 1.c1. ( r ) = r k + a1r k-1 + a2r k-2 + ... + ak = 0,
then Yn = rt is a solution to the homogeneous equation (3.9).
Proof:
Substituting Yn = rt into (3.9) give
n+k n+k-1 n+l n n k k-1
ri + a1ri + ... + ak-Iri + akri = ri (ri + a1ri + ... + ak-Iri + ak)
= rtf(ri) = rt.o = 0.
Hence, Yn = ri" is a solution to equation (3.9).
If the k roots of the characteristic equation f(r) = 0 are distinct, then the general
solution of (3.9) will be Yn = c1r1" + c2r2" + ... + ckrk"· In the case where ri has
multiplicity mi, i=1,2, ... ,j, where ml + m2 + ..... + mj = k, then the solution of Equation
(3.9) becomes
Yn=r1 "[ aii+ai2n+a13n2 + ...+ almlnm l-1] +r2" [ a21+a22n+a23n2 + ... ~ ~ ~ ' ~' ~
+a m2-1] "[ 2 mj-1] 2,m2n + ... + rj aj, 1 + aj,2n + aj3n + ... +aj,mP .
The solution Yn ofEquation (3.9) is usually written as YnH·
l__
41
A particular solution of the nonhomogeneous linear difference Equation (3.8) can
be derived from the relation y/ = f(Er1Rn. The general solution of Equation (3.8) then
becomes YnG = Yn" + YnP· When Rn has certain forms, we can use the trial solutions
given in Table 3.1 below. The unknown coefficients should be determined by
substituting the trial solutions into the difference equation. If a term of the trial solution
appeared in the homogeneous solution, then the entire trial solution corresponding to this
term should be multiplied by a positive integer power of n. This power should be just
large enough so that no term in the new trial solution will appear in the homogeneous
solution.
Example 3.3:
Solve Yn+2- 4Yn+1 + 4yn = 3(2)" + 5(4)"
Solution:
A homogeneous solution of this difference equation can be found by finding the
roots of the characteristic equation f(r) = r2 - 4r + 4 = 0. Since the roots are r1 = 2 and
r2 = 2, the homogeneous solution is y n H = c12" + c2n(2)". The particular solution for the
nonhomogeneous solution can be found by choosing the trial solution
y/ = A1n2(2)" + A2(4)". By substituting this y/ into the difference equation, we have
y/ = 3n2(2)"-3 + 5(4)"-1. The general solution ofthe difference equation is
Yn°= YnH + y/ = c12" + c2n(2)" + 3n2(2)"-3 + 5(4)"-1.
Example 3.4:
Solve numerically Yn+2 = -4Yn+1 + 5yn + h(4fn+1 + 2fn) and apply to
y I= y, y(O) = 1.
____...,__
_____....,___
Solution:
TABLE 3.1
A Trial Solution to Find a Particular Solution of Nonhomogeneous
Difference Equations
Terms in Tn Trial solution y/
~n A~"
sin an or cos an A cos an + B sin an
I
polynomial P(n) of degree m A on m + A In,m -1+ .. + A m
~"P(n) ~"(A0nm + A1nm-J+ .. +Am)
~n sm· an or ~n cos an ~"(A cos an+ B sin an)
-------
42
Calculating the polynomials p and cr, we have p(~) = ~2 + 4~ -5, cr(~) = 4~ + 2,
and p' (~) = 2~ + 4. The multistep method is consistent because p(1) = 1 + 4- 5 = 0 and
p' (1) = 2(1) + 4 = cr(1) = 4(1) +2 = 6. Since ao = -5, a1 = 4, a2 = 1, we have
ao + a1 + a2 = -5 + 4 + 1 = 0. By observation of formula (3.4), the condition is satisfied
for q=2 and q=3. It can be said then that the multistep method is of order 3. Since the
roots of characteristic p(~) = 0 are ~ 1 = -5, and ~2 = 1, then the multistep method is
unstable. If the computation satisfies the stability condition, then the method is said to
have convergence with order 3. The linear difference equation relation is
Yn+2 + 4(1-h)Yn+I- (5+2h)Yn = 0. As starting values, we can take Yo= 1, and y1 = eh. The
roots ofthe characteristic ofr2 + 4(1-h)r- (5+2h) are
r1 =2h-2+.J4h2 -6h+9 and r2 =2h-2-.J4h2 -6h+9.
l
60
40
20
Ul
~ 0
I
I
I
I
I
I
I
I
I
I • I
• II
I I
I ..
I I II
I I 1,1 I I
I
:>. 0.00 I
1.0tll
II
-20
-40
-60
tAxis
• Exact
• h=0.025
... h=0.05
.... h=0.1
l1
1
II
I\
I I
I 1
1
I I
I 1
1
I I
I 1
1
I I
I 1
1
I I
I 1
1 I.
Figure 3.1 Computational Results ofyn+2 = -4Yn+1+5yn+h(4fn+I+2fn)
Applied to y'=y, y(O)=l.
43
The general solution of the difference equation is Yn H = Ar1 n + Brt The computations
are done using h=0.025, 0.05, and 0.1 as shown in Figure 3.1. In this case, the smaller
the stepsize is chosen, the worse result is achieved. This occurred because as h ~ 0, y n H
will oscillate to infinity for r1 and r2 equal to 1 and -5, respectively. We can see that the
second part of y n H has a very strong growth in the solution for a small stepsize. If we
approximate the squareroot in r1 and r2 using the first two terms of Taylor's series, we
i J..
44
have r1 = 1 + h + O(h2), and r2 = -5 + O(h). Since r1 is an approximation of eh, the first
term of Yn H approximates the exact solution e1
. The second part of Yn H is often called a
parasitic solution.
3.2 Adams Methods
John Couch Adams was motivated by a problem given by F. Bashforth when
Adams proposed Adams-Bashforth methods. The Adams-Newton Backward Difference
formula is used as approximation functions to the integrands. If y is a function oft and t
is replaced by tk = t0 + kh, the backward difference of the Yk value is denoted by V'yk and
defined as Yk- Yk-1· It can be shown then that V'nyk = V'n-IYk - V'n-IYk-l· From the k known
values Yn-k+I> Yn-k+2' · · · ,Yn-1> Ym the k points (tu-kH,fn-k+l), (tn-k+2Jn-k+2), · · .,(tn,fn) can be
determined. The Newton Backward Difference formula interpolates the function f(t,y) of
y' =f(t,y) as
k-1 ( J p•(t) = p•(tn +sh) = L(-1)j -~ V'jfn
j=O J
(3.1 0)
(-sJ · s! t-t
where-. =(-1)J .
1
( • , s=-h", andt stst 1• Thesolutionofy'=f(t,y)at
J J.S-J)! n n+
t
n + 1
tn+l is given by y(t n+l ) = y(t n) + jf(t, y(t))dt, which in the numerical result is given by
t
n+l
Yn+l = Yn + fp* (t)dt
n
n
(3 .11)
By substituting (3 .1 0) into (3 .11) and using a new variable s such that the integration will
be from 0 to 1, Lambert [78] has shown that
k-1
- h"" • j Yn+l- Yn + ~y.V' f J n
(3.12)
j=O
where the coefficients Yi * satisfy
1
Y: = ( -1 )i J ( -~J ds
0 J
(3.13)
The multistep method (3.12) is called an explicit Adams-Bashforth method. Recurrence
relations for coefficients Yi* (see Appendix C) are given as
• 1. 1. 1.
Y. +-y. I +-y. 2 +-y. 3 +
I 2 1- 3 1- 4 1-
1 •
+ i + 1 y 0 = 1 , i = 0,1 ,2 ... (3.14).
The family of explicit Adams-Bashforth methods then can be written as
1 5 2 3 3
Yn+l = Yn +h(fn +lY'fn +
1
2 Y' fn +}8Y' fn + · · .). (3.15)
By considering the multistep method (3.15), special cases of (3.12) can be shown to be
k = 1 : Y n+ I = Y n + hf n.
k=2: Yn+l =Yn +~~fn -~fn-1).
,( 23 16 5 J k = 3: Yn+l = Yn +ll\__}2fn -}2fn-l + 12 fn-2 ·
,(55 59 37 9 )
k = 4: Y n+l = Y n + ll\__24 fn - 24 fn-1 + 24 fn-2 - 24 fn-3 .
These methods are very interesting because Yi*'s do not depend on k. They can
be determined directly from their recurrence relations (3.14). If we know the values of
46
y0*, y1*, y2*, y3*, we can create all k-step Adams-Bashforth methods with k=1, 2, 3, 4.
By observing Equation (3.12) or (3.15), we can see that from a k-step Adams-Bashforth
method, we can create a (k-1)-step Adams-Bashforth method just by throwing out the
last term of the series. We also can create a (k+ 1 )-step method by adding an extra term
after calculating coefficient 'Yk+ 1 * using their recurrence relations (3 .14 ).
Thus, we have already shown that the k-step Adams-Bashforth methods only
involve k terms (excluding y 0 ) of series(3 .12), and is given as
k-1
Y n+l = Y n + hL Y ~V'jf J n.
j=O
which has order k and error constant 'Yk *. By investigating the Adams-Bashforth
methods using test equation y ' = 'Ay, Hairer and Wanner [56] said that the Adams-
Bashforth methods are not appropriate for stiff problems.
Implicit Adams methods can be derived by interpolating the function f(t,y) of
y ' = f(t,y) over k points given in the explicit methods and one unknown point (tn+I,fn+I)·
In this case the approximation function off(t,y) becomes
k ( . -s+ 1 .
p(t)=p(t +sh)= '"'C-1)J )vJf n L..J J. n+l •
j=O
(3.16)
By substituting p*(t) in (3.11) with p(t), we have the implicit method
k
Yn+l = Yn + hL'"...J' 'V1 J. V'jfn +l" (3.17)
j=O
where the coefficients 'Yj satisfy
47
I
y; ~(-!);I (-s; l)ds. (3.18)
The multistep methods (3.17) are called implicit Adams-Moulton methods. Recurrence
relations for coefficients Yj (see Appendix C) are given as
1 1 1 1 {1 if i = 0
Yi +2Yi-l +)Yi-2 +4Yi-3 + · · · + i +1 Yo= 0 ifi = 1,2, .. (3.19).
The family of the implicit Adams-Moulton methods then can be written as
1 1 2 1 3
Yn+l = Yn +h(fn+l-2\lfn+I-UV' fn+l- 24\7 fn+l + · · .). (3.20)
By considering the multistep method (3.20), special cases ofEq. (3.17) can be obtained:
k = 1: Yn+l =Yn +{~fn+l +~fn).
k=2:
,(5 8 1 )
y n+l = y n + \12 fn+l +ufn -ufn-1 .
k=3:
,( 9 19 5 1 )
Y n+l = Y n + \24 fn+l + 24 fn- 24 fn-1 + 24 fn-2 .
k=4: = + 1 ( 251 f + 646 f - 264 f + 106 f - ___!2_ f )
Y n+l Y n \720 n+l 720 n 720 n-l 720 n-2 720 n-3 .
The same arguments as in the Adams-Bashforth methods also prevail in implicit
Adams-Moulton methods. But here, the k-step Adams-Moulton methods involve (k+l)
terms (excluding y 0
) of series (3 .17); That is, truncation is done after (k+ 1) terms
(excluding y 0 ) of series (3 .17). For k= 1, the backward Euler method y n+ 1 = y n + hfn+ 1 is
also called a 1-step Adams-Moulton method. It also can be shown that a k -step Adams-
Moulton has order (k+ 1) and error constant Yk+ 1. Shampine and Gordon [99] used these
...
..l
48
implicit Adams-Moulton methods using PECE (will be discussed later) m
developing ODEs packages to solve non-stiff differential equations. Byrne and
Hindmarsh used Adams methods to solve nonstiff differential equations in the EPSODE
[17].
The consistency, stability, and convergence of Adams-Bashforth and AdamsMoulton
methods can be examined using the stability concept proposed by Dahlquist.
Let us take an example for k=2. The explicit two-step Adams-Bashforth method is
Yn+1 = Yn + (h/2)(3fn - fn_ 1). By writing this formula in the form of Equation (3.3), we
have p*(r) = l - r and cr*(r) = (1/2)(3r - 1). Since p*(l) = 1 - 1 =0, and
p*'(1) = 2(1) - 1 = (112)(3.1-1) = cr*(1), then the explicit two-step Adams-Bashforth
method is consistent. The method is stable because the roots of the characteristic
equation p(r) are r1 = 0 and r2 = 1. By applying Equation (3.4) to this method, it can be
said that this method has order 2. The implicit two-step Adams-Moulton method is
Yn+1 = Yn + (h/12)(5fn+1 + 8fn - fn_ 1). As was done for the explicit two-step Adams-
Basforth method, we have p(r) = r 2 - r, and cr(r) = (1/12)(5r 2 + 8r - 1).
Since p(1) = 1 - 1 = 0, and p ' (1) = 2(1) - 1 = (1/12)[5(1)2+8(1)-1] = cr(1), then the
implicit two-step Adams-Moulton method is consistent. The method is stable because
the roots of the characteristic p(r) are r1 = 0 and r2 = 1. By applying Equation (3.4) to
this method, it can be said that this method has order 3
3.3 Backward Differentiation Formulae (BDF)
49
So far, we have discussed multistep methods based on difference equations and
integrations. We will now derive multistep methods which are based on numerical
differentiation of a given function. These methods are called the backward
differentiation formulas (BDF) and can be applied to solve stiff differential systems.
These implicit multistep methods are the first methods used in solving stiff differential
equations and become the most widely and predominantly used for all stiff computations
since Gear [44] proposed the computer codes in his book. The EPSODE created by
Byrne and Hindmarsh [17] also used these methods to handle stiff system problems. The
MEBDF created by Cash [20] to solve stiff problems is also an example of package that
implements BDF concepts.
The Newton Backward difference interpolations ofk points (tn+bYn+J), (tn,Yn), ... ,
(tn-k+I,Yn-k+J),(tn-k,Yn-k) are given as
k ( )
. -s+ 1 .
p(t)=p(tn+sh)=Pk(s)=~(-1/ j Y'Jyn+l" (3.21)
The values Yn-b Yn-k+b ... , Yn are known, and we would like to determine the value of
Yn+I· The methods are done by differentiating equation (3.21) and substituting the value
into y ' =f(t,y). The function of y is approximated by a polynomial p(t) of degree k.
Since the approximation ofy' = f(t,y) at tn+I is p' (tn+I) = fn+b we have (see Appendix
C)
1 d ~ ·(-s + 1) . --L,,(-1)J Y'J -
h ds j=O j Yn+l- fn+l'
which can be written as
..
k
Lbj'Vjyn+l = hfn+l"
j=O
By direct differentiating and substituting s= 1, the coefficients 8/ are obtained as
{
0, for j = 0
s: - 1 c. . .
u j - -: , 10f J ~
J
The multistep methods then can be written as
k
"L.".,. !. \7 j Yn +l -- hfn +l ·
j=l J
50
(3.22)
These multistep formulas are known as backward differentiation formulae (BDF
methods). Special cases ofBDF can be obtained:
k= 1 : Y n+ 1 = Y n + hfn+ 1 •
k=2: Yn+l = (4/3)Yn- (l/3)Yn-l + (2/3)hfn+l·
k=3: Yn+l = (18/11)Yn- (9/11)Yn-l + (2/11)Yn-2 + (6/11)hfn+l·
k=4: Yn+l = (48/25)Yn- (36/25)Yn-l + (16/25)Yn-2- (3/25)Yn-3 + (12/25)hfn+l·
k=5: Yn+l = (300/137)Yn- (300/137)Yn-l + (200/137)Yn-2- (75/137)Yn-3 +
( 12/13 7)y n-4 + ( 60/13 7)hfn+ I·
k=6: Yn+l = (360/147)Yn- (450/147)Yn-l + (400/147)Yn-2 -(225/147)Yn-3 +
(72/147)Yn-4- (10/147)Yn-5 + (60/147)hfn+l·
Since finding the roots of the polynomial p(r) of BDF methods (3.22) is not
simple, Definition 3.1 can not be used directly to investigate the stability of BDF
methods. By manipulating p(r), Hairer and Wanner [56] have shown that k-step BDF-methods
are stable only for k ~ 6.
51
3.4 Predictor-Corrector Methods
We have seen that the iterative procedures for the implicit multistep methods
require initial guess of the solution at each timestep. The simple way to handle this
problem is to use an explicit method to supply the required information. In this case, the
explicit method used to supply the information is called a predictor and the implicit
method used to correct the result is called a corrector. The combination of both
predictor and corrector is called a predictor-corrector (PC) method. The predictor-corrector
method consists of three processes: predicting (P) the next value using the
explicit method, evaluating (E) the derivative based on the latest value of y, and
correcting (C) the result using the implicit method. Normally, the algorithm of the
predictor-corrector methods are in the form of P(EC)m if they end with a correction or
P(ECtE if they end with an evaluation. The important guidelines [78] that should be
followed are:
A void choosing the order of predictor greater than that of the corrector,
because the local truncation error will still follow that of the corrector and
- Manage to choose the order of predictor equal to the order of the corrector, so
we can avoid unnecessary computations.
Suppose we use an implicit linear multistep method to solve the standard initial
value problemy' =f(t,y). An implicit linear multistep method can be written as
k-1 k-1
Yn+k + Lajyn+j = hbJ(tn+k'Yn+k)+hLb/n+j (3.23)
~0 ~0
and a linear multistep predictor as
52
k-1 k-1
Yn+k + La~yn+j = hLb~fn+j. (3.24)
j=O j=O
We are required to determine Yn+k from the formulas. Since we cannot solve Yn+k
directly, we can use the following procedure to calculate Yn+k:
P: Predict the value Yn+k(O) for Yn+k using (3.24).
E :Evaluate f(1n+k,Yn+k (0)) ·
C: Correct Yn+k(O) to obtain a new Yn+k(l) for Yn+k using (3.23).
(1) E :Evaluate f(1n+k•Yn+k ).
C: Correct Yn+k(l) to obtain a new Yn+k(2
) for Yn+k using (3.21)
The sequence of operations PECECE . . . to determine the value of Yn+k is given as
(0) (1) (2)
Yn+k ' Yn+k ' Yn+k ' · · · Lambert [78] proposed a single mode P(ECtE1-t instead of
two modes P(ECtE or P(EC)m as follows:
P(EC)mE1-t:
k-1 k-1
P: [OJ +"a~ [m] = h "b•f[m~t]
Yn+k L. JYn+J L. J n+J
j=O j=O
[v] [v]
fn+k = f(tn+k,Yn+k) l
I
(ECt: ~ v = 0,1, ... ,m-1
k-1 . k-1 1
[v+1] +"a. [m] = hb f[v] + h "b.f[m~t]J
Yn+k L. JYn+J k n+k L. J n+J
j=O j=O
E(l-t). f[m] = fi(t Y[m] ) ift = 0
· n+k n+k' n+k '
53
where t = 0 or 1. He also has shown that if p and p* are the order of corrector and
predictor respectively, then
1. if p*~p (or if p* < p and m > p-p*), the PC method and the corrector have the
same order and the same principal local truncation error,
2. if p*<p and m = p-p*, the PC method and the corrector have the same order
but different principal local truncation error, and
3. ifp* < p and m~p-p*-1, the order of the PC method is p*+m (<p).
Based on Lambert [78], most of modem predictor-corrector codes for non-stiff
problems use Adams-Moulton methods as correctors and Adams-Bashforth methods as
predictors. These methods are sometimes called ABM methods. Considering the
principal local truncation error given above, it is important then to choose the predictors
and correctors with the same order. To satisfy this condition, the ABM methods should
be:
k·l
Yn+l =Yn +hLy.j'Vjfn , p*=k.
j=O
k-1
Y n+l = Y n + hLy j 'Vjfn+l ' P = k
j=O
3.5 Extrapolation Methods for Solving ODEs
The meaning of extrapolation in the sense of integration methods is the process of
deriving an improved estimate for the value of a definite integral based on two or more
applications of a formula using different stepsizes of integration. The following example
...l
is an illustration of the meaning of extrapolation. Assume Yn(h) is the approximate
solution of y ' = f(t,y) which is integrated using stepsize h, then Yn can be written as
Yn(h) = y(tn) + hP E(tn) + O(hp+1
), where y(tn) is the true solution, E(t) is called the
magnified error function. We would like to manipulate h so that we have a more accurate
estimate to the true value ofy(tn). Let us take a new stepsize Ah (0 <A <1), and apply to
the formula ofy0 • Then, we have Yn(Ah) = y(tn) + (Ah)P E(tn) + O(hp+1
), By multiplying
Yn(h) with -AP and then adding to Yn(Ah), we have
Yn(Ah)- APYn(h) = (1 - Ap) Y(1o) + O(hp+l),
and also
Yn(Ah) + [Y0(Ah)- Yn(h)]/[(1/A)P- 1] = y(tn) + O(hp+1
).
Ifwe choose Yn(Ah) + [y0 (Ah)- Yn(h)]/[(1/A)P- 1] as the approximate value ofy(tn), the
error becomes smaller. This approximate value is called extrapolation or, sometimes
Richardson's extrapolation. Assume that a certain numerical formula has a trurJcation
error proportional to h3
. Suppose a computation of this formula with h = 0.10 gives an
output 4.5800 and with h = 0.05 gives an output 4.4350. What is the best approximation
to the true value based upon this information? In this problem A= 1/5 and p = 3. Using
the formula given above, the approximate solution is
y(0.05) + [y(0.05)-y(0.10)]/[(1/5)3-1] = 4.4350 + [4.4350- 4.5800]/[124]
= 4.4288.
In most of implementation, usually the factor A is chosen 1/2, so the Richardson's
extrapolation at point to is
55
h
(
h) Y n(-)- Y n(h)
y(t ) = y - + 2 + O(hP+1)
n n 2 2P -1
(3.25)
We have seen that the extrapolation depends on the time step h. This procedure can be
done repeatedly until h ~ 0 so that the approximate value tends to the exact solution.
Deuflhard [32,33] has shown some extrapolation methods including control order and
step size.
Let y(t) be the true solution of y ' = f(t,y), y(t0) = y0, and y(t,h) be an
approximate solution satisfied by choosing the stepsize h to the appropriate numerical
method. Assume y(t,h) can be written as an asymptotic error expansion in powers ofh,
r r r r r
y(t,h) = y(t) + a 1h 1 +a2h 2 + a 3h 3 + ... +amh m + am+lh m+ 1 + ...
where 0 < r1 < r2 < ... < rm, a~> a2, a3, ... are independent of h and determined by
evaluating y(t,h) over stepsize hi> i=O, 1, 2, . . . One approach can be obtained by taking
ri = ir and the stepsize sequence hi such. t hat ho > h1 > h2 >. . . Substituting values ri, we
have
y(t,h) = y(t)+alhr +a2h2r +a3h3r + ... +amhmr + ...
If the superscript k is related to the approximate value with stepsize hk and the subscript
m is related to the number of times the linear combination has been applied to eliminate
a~> a2, ... , ~' and Y m (k) is defined as
hr y(k+l) _ hr y(k)
y(k) = k m-1 k+m m-1
m h~-hk+m
(3.26)
then we will have
y::> = y(t)+(-l)mh~h~+lh~+2 ... h~+m(am+l +O(h~)).
56
Form ~ oo, Y m(k) will close to y(t). The representation of this extrapolation method is
shown in Table 3.2. The Table is expressed in a graphical table as
y (k)
m-1
~
y (k+l)
m-1
Yo tu, = y(t,ho)
Yot'' = y(t,hi)
y~ =y(t,h2)
YotJJ = y(t,hJ)
Yot4J = y(t,h4)
y (k)
m
TABLE3.2
Table of Extrapolation
y
1
tUJ
ylt~J y 2lUJ
y
1
tLJ YtJ y
3
tUJ
YtJ y2t.l) y3t1J
Equation (3.26) can be written as
y 4lVJ
y(k) = y(k+l) h~+m (Y~~/) - hkr y(k) )
m m-1 + +m m 1
hr - hr
k k+m
For h=ho:A} with 0 <A.< 1, we have
(3.27)
y(k) = y(k+l) bmr (Y(k+l) - y(k) )
m m-1 + m-l m I
1- bmr
For b= 1/2, we have
y<k+1) - y(k)
y(k+1) + m-1 m-1 ifr = 2
m-1 2m 1 '
y(k)-
m -
y(k+1) y(k+1) - y(k)
m-1 + m-1 4m -1 m-1 ' ifr 4
Equation (3 .28) is sometimes called Romberg integration.
57
(3.28)
Let us take an example y ' = -y, y(O) = 1. We will interpolate this problem until
t = 1 using Euler extrapolation with r = 2 and the number of points in the interval chosen
to be m = 10. After that we approximate the problem using Euler extrapolation with
m = 5, Euler method and ABM method of order 4 with h = 0.25, and integrate until
t = 2.5. For the first Euler extrapolation, we choose hi= (1-0)*2-i, i=0,1, ... ,9. Using
Euler method, we have a relation Yn+l = Yn + hf(tn,yn) = (1-h)Yn· Substituting y0, we have
a relation y n = ( 1-h)n. To reduce the use of storage, we can use variable y with dimension
as much as m. A procedure of the Romberg extrapolation is represented as:
Y
-y (0) o- o
Y -y (I)
1- 0
Y
-y (2)
2- o
Y1= Y1 + [yi-Yo]/[2
1
-1]
Y2= Y2 + [yz-yJl/[2
1
-1] Y2= Y2 + [yz-yJl/[i-1]
Y(m-I)=Yo (m-l) Y(m-1) = Y(m-I)+[Y(m-lrY<m-2)]/[2
1
-1] Y(m-1) = Y(m-1) + [Y(m-lrY<m-2)]/[2
2
-1]
,1
:j
An algorithm of Romberg integration using Euler method is given as:
1. h0 = (b-a) I* a and bare starting and final times respectively *I
2. Compute y(i), i = 0, 1, 2, ... ,m-1
for i = 0 to m-1 do
end do
h = ho*2-i
n=2i
y(i) = (1-h)n I* With stepsize h, y(l) is achieved after
n=2i steps *I
3. Extrapolate the values
for i = 1 to m-2 do
k=m-1
end do
for j=i+ 1 to m-1 do
end do
y(k) = y(k) + [y(k)- y(k-1)]1[2i- 1]
k =k-1
58
The computation table of the algorithm for Romberg integration using Euler method
given above for the case ofy' = -y, y(O) = 1 is given in Table 3.3.
The table is derived as follows: Let us take m = 4 and i =5.
y (i) = y (i) + [Y (i) _ y (i-1)] I [2m-1]
m m-1 m-1 m-1
=Y/5> + [Y3<5>- Y3<4>] I [24-1]
59
= 0.3678811228 + [0.3678811228-0.367920598]/15 = 0.367878603.
TABLE3.3
Romberg Representation ofy' = -y, y(0)=1
1 h m Yi,m=O Yb m=1 Yb m=2 Yi,m=3
0 ho2v 2v 0.0
1 ho2-' 2' 0.25 0.500000000
2 h0T'' 2" 0.316406250 0.382812500 0.343750000
3 h0T' 2j 0.343608916 0.370811582 0.366811275 0.370105743
4 ho2_.. 2'+ 0.356074130 0.368539345 0.367781933 0.367920598
5 h0T'' 2" 0.362055289 0.368036448 0.367868816 0.367881228
6 hoTo 20 0.364986524 0.367917759 0.367878196 0.367879536
7 hoT 2' 0.366437716 0.367888908 0.367879290 0.367879447
8 ho2-ll 21:1 0.367159755 0.367881794 0.367879423 0.367879442
9 ho2-~ 2~ 0.367519891 0.367880028 0.367879439 0.367879441
I Yb m=4 . [_ Yi> m=5 I y· m=6 I' I Y· m=7 I'
60
0.367774922
0.367878603 0.367881947
0.367879424 0.367879441 0.367879410
0.367879441 0.367879441 0.367879441 0.367879441
0.367879441 0.367879441 0.367879441 0.367879441
0.367879441 0.367879441 0.367879441 0.367879441
For the second computation of Euler extrapolation, we choose m = 5, with
stepsize for all methods h = 0.25. The computation results is given in Table 3.4. It can
be shown that we can always control the accuracy of results produced by extrapolation
methods, just by controlling the stepsizes over the entire range of computations.
Graphical representations of these computations compared to those resulted by Euler
method and Adams method of order 4 are shown in Figure 3.2. The error vs timet of
this extrapolation method tends to be accumulated as shown in Figure 3 .2.
t
0.0
0.25
TABLE3.4
The Computation Results ofy' = -y, y(0)=1
Using Extrapolation, Euler, and ABM-4 Methods
Exact Extrapolation Euler ABM-4
1. 000000000 1. 000000000 1. 000000000 1.000000000
0.778800783 0.778800747 0.750000000 0. 778808594
l_
0.50
0.75
1.00
1.25
1.50
1.75
2.00
2.25
2.50
rn ·x
1.00
co 0.50
>.
0.25
0.606530660 0.606528620
0.472366553 0.472345756
0.367879441 0.367774922
0.286504797 0.286148404
0.223130160 0.222179660
0.173 773943 0.171635077
0.135335283 0.131086662
0.105399225 0.097729534
0.082084999 0.069250702
Exact
Extrapolation
Euler
.. ABM-4
0.06
0.05 -
0.562500000
0.421875000
0.316406250
0.237304688
0.177978516
0.133483887
0.100112915
0.075084686
0.056313515
0.606542826
0.472380765
0.385532323
0.299920011
0.233222895
0.181836888
0.141601456
0.110261898
0.085867179
Extrapolation
Euler
ABM-4
/~ . ~
0.04 ,_ I "'"' we 0 .03·- r ~-
0.02
0.01 /
0.00--- ------ ------- • • 0•.5 • c~~~~ 0~ ~5 1~ 1~ 2~ 2~ 1.0 1.-5 2 .O 2.5
taxis taxis
Figure 3.2 Computation Results ofy' = -y, y(0)=1 Using
Extrapolation, Euler, and ABM-4 Methods
61
..L
CHAPTER IV
STIFF ORDINARY DIFFERENTIAL EQUATION
Dealing with system of differential equations will sometimes lead us into
stiffness problem. This possibility occurs especially when the solution contains a fast
mode which decreases very fast in response to initial conditions or a changing input and
a slow mode which needs longer computing time for its transient to decay. This situation
commonly occurs when standard numerical methods are implemented to find the
approximate solutions of a differential equation having the exact solution containing a
term of the form e"\ where A is a complex number with negative real part. As
mentioned before, this term will tend to zero when t increases to infinity, but
unfortunately its approximation generally will not always have the same property, except
when the stepsize is chosen such that the method used is still stable. The problem will
become serious when the exact solution has a steady-state term that does not grow
significantly with t, and a transient term that decays fast to zero. In this case, the
numerical methods can approximate the steady state solution, but they need attention
when they are applied to approximate the decaying transient term that can dominate the
computation, and can produce incorrect results. A standard numerical method can be
used to solve stiff problems, but the cost will be very high since we have to maintain an
62
63
extremely small stepsize for the entire of integration [13,105]. The two following
problems will help us to understand the stiffness concepts in ordinary differential
equations.
Problems 4.1:
Consider the following set of equations
( y;J ( 0 1 )(y') (y,(O)) ( 1)
y~ = -100 -101 y2 ' Y2(0) = -1
Discussion:
The eigenvalues of the matrix
A= (-~00 -1
1
01)
are A = -1 and A= -100. These eigenvalues and the initial conditions give solution
(
y1(t)) ( e-1 J
y2(t) - -e-1
From this exact solution, we cannot see any problem at all, but when the problem is
solved numerically, the unstable computation arises unless a small stepsize is used. The
corresponding solution resulting from the eigenvalues are e-t and e-Ioot as shown in Figure
4.1. The initial condition has thrown the e-IOOt term out from the exact solution. By
considering the eigenvalues, we know that the dependent variables contain both fast and
s1 o w components. I n t hI.S case, e- lOOt de cays f:a st to zero, wh "1l e e- t de creases to zero very
slowly. To maintain the stability of a simple Euler method, the stepsize should be
chosen such that 11 +Aminhl ~ 1. Therefore, the stepsize of the computation should be
maintained within the intervall1-1 OOhl ~ 1 or 0 ~ h ~ 2/100 = 0.020 .
•
64
1
"l
I
1
~
I
i
I
tJ)
~ I \ ' e·t
>-
~
'\
\
\
\ 9 -1oot
0 '
0 2
tAxis
Figure 4.1 Graphical Solutions of e-t and e-Ioot
Problem 4.2:
Consider the following set of equations
(y·;) = ( 1 2 )(Y1) , (y1(0)) = ( 1)
y2 -4 -5 y2 y2(0) -1
Discussion:
The eigenvalues of the matrix
A~(~ :5)
65
are A= -1 and A= -3. The corresponding solutions resulting from the eigenvalues are e-t
and e-Jt respectively. These eigenvalues and the initial conditions give solution
(
y1 (t)) ( e-1
)
y2(t) - -e-1
Using Euler method, the absolute stability can be maintained if 11 + Aminhl ~ 1. Thus, the
stepsize of computation should be maintained in the intervall1-3hl ~ 1 or 0 ~ h ~ 2/3.
Both problems have the same exact solutions, but they behave very differently
when they are solved numerically. In this case, we need to apply different stepsizes to
compute the two problems using Euler method. The cost of Problem 4.1 is far more than
that of Problem 4.2 when they are solved numerically. The cost for integrating Problem
4.1 from t0 = 0 to tr = 50.0 is at least 2500 steps (for h = 0.02), while the cost for
integrating Problem 4.2 from t0 = 0 to tr = 50.0 is at least 75 steps (for h=2/3). The
phenomenon illustrated in these two problems is known as stiffness. Problem 4.1 is stiff,
while Problem 4.2 is non-stiff. Since both problems have the same exact solutions, the
phenomenon should be a property of the differential system itself.
Even though the exact solution for eigenvalue A = -100 of Problem 4.1
contributes practically nothing to the solution of y1 and Y2> the criterion of absolute
stability forces us to use an extremely small value of h over the entire range of
integration. Therefore, it can be said that using the same numerical method, the
computation time needed to solve a highly stiff differential system is more expensive
than that of non-stiff one. The main objective of a good ODE solver is to minimize the
computing cost with subject to the required tolerance and the stability of computation. In
...
66
order to gain this objective, an ODE solver is usually equipped with tools for: detection
stiffness, determining a starting stepsize, and controlling stepsizes and formulae, by
considering the required tolerance.
4.1 Stiff Differential Problems
The study of stiff systems has become popular in the last twenty years because
they arise in many applications such as: radio technology, chemical kinetics, reactor and
radiation physics, hydrodynamics, process dynamic and control, electrical circuits,
diffusion, beam, and lasers. Several methods have been proposed and analyzed
carefully, but only few of the codes have been developed until recently. In this thesis,
the investigation will be done based on stiff differential equation packages such as
MEBDF, LSODE, EPSODE, and VODE, using the test cases:
1. Kidney problems introduced by Scott and Watts in the form [17, 21, 79]
Y2Y1' = ay1(Y3- YI),
Y2' = -a(y3- YI)
Y4Y3' = b- c(y3- Ys)- ay3(y3- YI)
Y4' = a(y3- YI)
YI(O) = 1
Y2(0) = 1
Y3(0) = 1
Y4(0) = -10
Ys' = -c(ys- Y3)/d Ys(O) =A
where a= 100, b = 0.9, c = 1000, d = 10,0 :$; t :$; 1, with initial conditions
Problem A
Gl 0.9902688359
G2 0.9902834990
03 0.9925211341
04 1.0304879856
05 0.99
06 0.9
07 0
2. An autocatalytic reaction pathway proposed by Robertson [21, 56]
Yl' = -0.04 Yl + 104 Y2Y3,
Y2' = 0.04 Yl- 104 Y2Y3- 3.107 y/,
Y
,_
3 -
for 0 ~ t ~ 4. 108
.
7 2 3.10 Y2 ,
3. Problem D4 of Enright et al. [39, 103]
Yl' =- 0.013yl- 1000YIY3
y2' =- 2500y2 y3
y3' = 0.013y1 - 1000 y1y3 - 2500y2y3
where 0 ~ t ~ 50.
4. Problem proposed by Gupta and Wallace [53]
Y 1' = vy 1 - wy 2 + (-v + w + 1 )e 1
Y2' = wyl + vy2 + (-v -w + 1)e1
The exact solution is
y vt ( t 1 = c1 e cos wt + c2) + e
y vt • ( t 2 = c1 e sm wt + c2) + e
YI(O) = 1
Y2(0) = 0
yJCO) = 0
Y1(0) = 1
Y2(0) = 1
Y3(0) = 0
YI(O) = 1
Y2(0) = 1
67
68
where v = -80, w =8, 0$ t $ 10.0.
4.2 Stiffness Concepts
There are two approaches that have been taken in order to define stiffness. The
first approach is to define stiffness based on various aspects of phenomena of stiffness.
The second approach is to define stiffness in mathematical terms. Even though several
definitions relating to mathematical terms of stiffness have been proposed, but they are
all still in debate until now. Stiffness can not be defined precisely in mathematical terms
[1, 62, 78], even for the class of linear constant coefficient systems. The following
differential equation concepts are useful in order to understand the stiffness concepts.
The homogeneous solutions of ODE problems are also called transient solutions,
and the particular solutions are called steady-state solutions. For the cases of linear
systems, the homogeneous solutions correspond to the characteristic roots A.i of their
matrices. If Re(A.i) < 0 and IRe(A.i)l is large, the corresponding homogeneous solution
will decay fast as t increases and is thus called a fast transient, but if the IRe(A.i)l is small,
the corresponding homogeneous solution will decay slowly and is called a slow transient.
Suppose ~and~ in {1-i ,i = 1, .. ,n} are defined by IRe(~)l $jRe(A.)j $1Re(~)l, Vi. We
will have trouble when IRe(~)l is very large and jRe(~)j is very small, because we are
forced to integrate over the whole range of integration with a very small stepsize.
69
The following discussion of stiffness concepts will be based on Lambert's
explanations [78]. If we take the ratio IRe(~)I!IRe(~)l as the stiffness ratio, the statement
as a candidate for the definition of stiffness is given as follows:
Statement 1:
A linear constant coefficient system is stiff if all its eigenvalues have negative
real part and the stiffness ratio is large (> > 1 ).
This statement has been used by Lambert [1973] as a definition for stiffness, but now he
realizes that the definition is not correct any longer. He used the following three systems
to show that the definition is not satisfactory with the phenomenon.
System 4.1:
(y~J (-2 1 xy1J ( 2sint J ly~ = 998 -999 y 2 + 999( cost- sin t)
The general solution of this system is given as
(y 1 (t)J = c e -t(1J + c e -woot( 1 J +(sintJ
y2 (t) 1 1 2 -998 cost
The stepsize interval of stability of this system using an Euler method is
0$ h $2/1000.
System 4.2:
( y~ l ( -2 1 xy1J ( 2 sin t J
lyJ= -1.999 0.999 y2 + 0.999(sint-cost)
The general solution of this system is given as
(
yl(t)J=ce-t(1J+ce-o.ooit( 1 ]+(sintJ
y 2 (t) 1 1 2 1.999 cost
70
The stepsize interval of stability of this system using an Euler method is 0~ h ~2.
System 4.3:
( y~ l (-0.002 0.001 xy1J ( 0.002 sin(O.OOlt) J
lyJ = 0.998 -0.999 y2 + 0.999[cos(O.OOlt)- sin(O.OOlt)
The general solution of this system is given as
(
Y1 (t)J = -t( 1 J -o.ooJt(1J (sin(0.001t)J
y
2
(t) c1e -998 + c2e 1 + cos(O.OOlt)
The stepsize interval stability of this system using an Euler method is 0~ h ~2.
It is easy to show that the three systems have the same stiffness ratios: 1000.
Even though they have the same stiffness ratios, System 4.1 and System 4.2 have
different stepsizes in order to maintain the absolute stability of computations. No matter
how we choose an initial condition for System 4.1, the numerical method needs a very
small stepsize for integration. The effect of stiffness is that we have to continue using
that small stepsize long after the fast transient solution has no influence. In this case,
System 4.1 is called a genuine stiff [78]. Suppose we equip System 4.2 with initial
conditions: y(O) = [2,3.999]T causing c1 = c2 = 1. It is obvious that if we want to reach
steady state solutions, the transient solutions should vanish. For System 4.2, this is only
poss1" bl e 1"f t h e so 1u t1. 0n e -O ·O O!t 1• s c1 o se to zero. In t h"1 s case, t h e tota1 computati.O n step of
System 4.2 is comparable to that of System 4.1, because we still need to compute e-o.ooit
until it has no effect on the solution of System 4.2. In this situation, the definition of
stiffness in Statement 1 seems to be consistent. The conclusion will be different when
the initial conditions for System 4.2 is changed to be: y(O) = [2,3]T causing c1 = 2, and
...
71
c2=0. Using these initial conditions, the slow transient term disappears from the exact
solution so there is no need to integrate a long way until the steady state condition is
reached. From these two initial conditions, we see that the definition of Statement 1 is
not consistent any longer, because it depends also on initial conditions imposed by a
particular problem. In this situation, System 4.2 is not genuinely stiff, but we can called
them pseudo-stiff. By doing the same evaluation, it can be said that System 4.1 is also
pseudo-stiff. System 4.3 can be changed to System 4.1 by making the transformation
t = lOOO't. In this case, System 4.3 is stiff in exactly the same sense as System 4.1.
Statement 2:
Stiffness occurs when stability is more of a constraint than accuracy.
This statement also fails, because the stability also depends on the accumulation of
errors. By examining the local error at the very first step, it can be said that a stiff
problem has a substantially higher local error than that of a non-stiff problem.
Statement 3:
Stiffness occurs when some components of the solution decay much more rapidly
than others.
This statement also fails, because it is merely based on comparison of the rate of change
of the fastest transient solution with that of the steady-state solution. From the discussion
of Problem 4.1, this statement fails totally.
Definition 4.1:
If a numerical method with a finite region of absolute stability, applied to a
system with any initial conditions, is forced to use in a certain interval of
L
72
integration a stepsize which is excessively small in relation to the smoothness of
the exact solution in that interval, then the system is said to be stiff in that
interval.
This definition can distinguish between the genuine stiff and the pseudo-stiff systems,
and it offers the idea that stiffness can be considered to change over the interval of
integration. To get a good understanding of this definition, we have to do some
computations, otherwise it will be difficult to understand.
Statement 4:
A system is said to be stiff in a given interval of t if in that interval the
neighboring solution curves approach the solution curve at a rate which is very
large in comparison with the rate at which the solution varies in that interval.
The same case with Definition 4.1 occurs in Statement 4.
4.3 Stability Concepts of Numerical Stiffness Methods
Several definitions of stability related to stiff systems are given in order to satisfy
the needs of smaller stepsize more for stability than for accuracy.
Definition 4.2:
A method is said to be A-stable in the sense of Dahlquist if when it is applied to
the test problemy' = 'Ay, the solutions converge to 0 as n ~ oo or the region of
stability is a superset of { 'Ah I Re('Ah) < 0}. A minimum region of this stability is
shown in Figure 4.2.
lm(A.h)
A-Stable
0 Re(A.h)
Figure 4.2 A representation of a Minimum
Region of A-stability
73
It can be shown that a trapezoidal formula, a second-order Adams-Moulton
method Yn+I = Yn +(h/2(Yn+J 1 + Yn1
), is A-stable. By substituting the test function y 1 = 'Ay
into the trapezoidal formula, we have Yn = [(2+/...h)/(2-'Ah)t y0. The region of stability is
reached if 12 + 'Ahl/12 - 'Ahl ~ 1. Here, the region is the set of points z such that the distance
between z and the point z = 2 is greater than the distance between z and the point z = -2.
This region is { z I Re(z) ~ 0}. That is, the trapezoidal formula is A-stable. The backward
(implicit) Euler method Yn+I = Yn + hYn+J 1 is also A-stable (see Chapter 2.1).
Definition 4.3:
A method is said to be A( a)-stable, a E (O,n/2) in the sense ofWidlund if when it
is applied to the test problemy 1 = 'Ay, the solutions converge to 0 as n ~ oo with
h fixed for all !arg('Ah)l < a, !'AI "# 0 or the region of stability is a superset of
{'Ah I -a < n - arg('Ah) <a} ; it is said to be A(O)-stable if it is A(a)-stable for
some a E (O,n/2). Figure 4.3 illustrates a minimum region of this stability.
Definition 4.4:
lm(A.h)
Re(A.h)
Figure 4.3 A Representation of a Minimum
Region of A( a)-Stability
74
A method is said to be A0-stable if the region of stability is a superset of
{A.h I Re(A.h) < 0, Im(A.h)=O}. Figure 4.4 illustrates a minimum region of this
definition.
Im(A.h)
+---------....j-;:----.Re(A.h)
0
Figure 4.4 A Representation of a Minimum Region of A0-Stability
I
~~
,,'
75
Definition 4.5:
A method is said to be stiffly stable ifR1uR2 is a subset ofthe region of stability,
where R1 = {.Ah I Re(.Ah) <-a} and R2 = {.Ah I -a~ Re(.Ah) < 0, -c ~ Im(.Ah) ~ c},
and a and c are positive real numbers. Figure 4.5 illustrates this definition.
Definition 4.6
Im(A.h)
I
RI
lc
R2
-a ----to Re(A.h)
-c
Figure 4.5 A Representation of a Minimum
Region of a Stiffly-Stable Method
A one-step method is said to be L-stable if it is A-stable and, in addition, when
applied to the scalar test equation y' = A.y, A. a complex constant with Re(A.) < 0,
it yields Yn+I = E(.Ah)y"' where IE(A.h)l ~ 0 as I.Ahl ~ -oo.
From the definitions, the hierarchy of stability is drawing as follows: L-stable ~ A-stable~
Stiffly-stable~ A(a)-stable ~ A(O)-stable ~ A0-stable.
4.4 A Modified Euler Method for Solving ODEs
76
Since most of the numerical methods used to solve stiff differential problems are
implicit methods, it is almost impossible to have a cheap low-accuracy code for stiff
ODEs. The meaning of cheap here is that the methods are easy to program and the
computing time for solving the problems is not too expensive. It is true that if one
persistently solves stiff problems, one should use one of the better codes which are now
accessible. But if one needs to compute stiff problems where the accuracy is not taken
into consideration especially for a quick low-accuracy solution, Lambert [77] proposed a
method to satisfy that requirement.
The idea that Lambert suggested is to apply an automatically determined
sequence of variable stepsizes with Euler's rule. Unfortunately, there is no guarantee
i
·~
that this method can be successfully used in order to satisfy the meaning of cheap.
Lambert claims that this method is stable and the solution is not unimportant, but
sometimes the computation is trapped into unacceptably slow progress. Even though
Euler's rule using the automatic sequence of stepsizes seems not to be a likely method
for solving stiff systems, this method can provide a usable solution with an acceptable
computing time for some stiff problems.
.;
Here, this method is introduced in order to have a method for solving stiff
problems that is easy to program and gives relatively lower computing time. The test
problems that will be used are:
1. Y1' = -YI - 0.5y2- 0.5 Y3 YJ(O) = -1
y2' =- 0.5y1- 1000.75 Y2 + 999.25 Y3 Y2(0) = 1
Y3' = -0.5y1 + 999.25y2- 1000.75y3 Y3(0) = 3;
2.
A.1 = -2000, A.2 = -2, A.3 = -1/2
with the exact solution
y 1 (t) = eA.2t - 2eA.3t
y 2 (t) =-eA. It + eA.2t + eA.3t
y 3 (t) = eA.It + eA.2t + eA.3t
y11 = 0.01- [1+(y1 +1000)( Y1 + 1)](0.01+ Y1 +y2)
Y21 = 0.01 - (1 + y/)(0.01 + Y1+ Y2)
77
Y1(0) = 0
Y2(0) = 0
The derivation of the method for solving y 1 = f(t,y), y(t0) = Yo is done by
introducing a parameter 8, such that the method called the 8-method, and is written as
Yn+1 = Yn + h[(l-8)fn+1 + 8fnl· The parameter 8 then should be determined such that the
method is exact when applied to the test problemy'= ll-Y· By modifying Euler's rule the
method becomes Yn+1 = Yn + hvf"' where v = 1+0(h) to guarantee the consistency of the
method, and our task then is to find v such that the modified Euler method holds exactly
when applied to the test problem y 1 = ll-Y· The exact solution of the test problem is
y = y0eJ.l1
• After substituting y 1 = lJ.Y into the modified Euler's rule, we have
v = [(Yn+1/Yn)-1]1hlJ. = [ehJ.l_1]1hlJ., where Yn+1 = y0e(n+1)hJ.l and Yn = e"hJ.l. Since the Taylor's
series of ehJ.l = 1 + hlJ. + O(h), we have v = 1 + O(h).
Consider the initial value problemy 1 = f(t,y), y(t0) = y0, y,f E Rm and the test
problemy 1 = Ay, y(t0) = y0, where A is a constant m x m matrix with real eigenvalues "-t
, t = 1,2, ... ,m satisfying A.1 ::;; A.2 ::;; A.1 ::;; .... "-m-1 ::;; Am::;; 0. By adopting the idea given
above, for the ODE system where An = Bf/By!"' we have
,J
'1
78
Yn+l = Yn +hv n f n 1
v, +'"· -l]ihJ.l, f. (4.1)
Jln =(f:Anfn)l((fJJ
On the irregular-spaced point set { tn I tn+I = tn + hn , hn = hvn }, the equation (4.1)
becomes
Yn+l = Yn + h n f n 1
h = [e efln - 1] I ~. n lln
1
(4.2)
lln = (f: Anfn) I C(fn)J
Since Equation (4.2) may have a large step, which would not be good for the accuracy,
Equation ( 4.2) is modified as
Yn+l = Yn + h n f n
hn = min(h0 , h:)
h: = [ e Sfln - 1] I lln
1
~.
I
lln = C( AJJ I (f)n)j
(4.3)
The value e should be chosen positive so that hn * is always positive. Since nothing is lost
if e is chosen very large, the modified Euler method with the automatically determined
sequence of stepsizes can be written as
,,'
l
!i
il
I
'
79
Yn +l =yn+ hnfn lI
h, = min(h,,h:) f (4.4)
h: = (f)n)ll( Anfnlj
4.5 An Explicit Exponential Method for Solving ODEs
Even though some authors define stiffness equations as problems that cannot
be solved by explicit methods [56, 62], there are others that have still tried to solve
stiffness problems using explicit methods. Among them are Treanor [79], Lawson [79],
Osborn [79], Lambert [77], and Ashour and Hanna [3]. The effort for solving stiff
problems using explicit methods is still worthwhile because implicit methods need space
and the cost of computing time per step is relatively higher than that of explicit methods.
Treanor [79] modified the fourth-order Runge-Kutta method for solving stiff problems,
Lawson [79] proposed a transformation of the stiff system into a nonstiff system, and
then applied a Runge-Kutta method, Osborn [79] used Pade fractions to approximate ehA.
, Lambert [77] modified Euler method (Chapter 4.4) and Ashour and Hanna [3] used an
exponential method by applying Watson's lemma of asymptotic expansion of Laplace
integrals. The following explicit exponential method was a method that was proposed by
Ashour and Hanna to solve stiff ODEs.
The explicit exponential method is developed based on the idea of detecting any
exponential decaying variables at each step of the computation and to approximate them
such that they are not exponential decaying any longer. To achieve this idea, Ashour and
80
Hanna [3] applied Watson's lemma. We know that integration of y ' = f(t,y) from t to
t+h is written as
t+h
y(t +h)= y(t) + Jf(s,y(s))d . (4.5)
Computations will decay exponentially if the local Jacobian fy(to,y0) is negative and
lfy(t0,y0)1 >> 1. If this situation occurs, the solutions are approximated with functions
such that computations tend not to decay exponentially any longer. The following
expansion in Laplace integrals will be used to implement the idea of choosing solution
methods aforesaid,
A
I= Je-15F(s)ds, t > 0 (4.6)
0
where F(s) is of exponential order. It means that function F is of exponential order ~ if
there exist positive real numbers T, c, and ~ such that IF(s)l < ce~5 , s > T. Based on
Watson's lemma, the integration (4.6) can be approximated by expanding F(s) using a
Maclaurin series and integrating it. In the case where the upper limit of A is too small
relative tot, we use an expansion that is asymptotic either for t~oo with A fixed, or for
A~O with t fixed. Using a mathematical manipulation involving the Maclaurin series,
we have
F(s) = e·msF(s)ems
= e-ms [F(O)e0 + s{F'(O)e0 + mF(O)e0
} + O(s2
)]
= e-ms[F(O) + s{F'(O) + mF(O)} + O(s2
)].
l
j
'
81
Let us take {F'(O) + mF(O)} = 0 such that the O(s) term in a Maclaurin expansion is zero.
By substituting the latter results, from ( 4.6) we have
I= }e-<t+m}sF(O)ds = -F(O) [e·(t+m}A- eo]= F(O) [1- e·(t+m}A].
t+m t+m
(4.7)
0
By taking s = t + u, Equation ( 4.5) can be transformed into
u=h
y(t +h)= y(t) + Je-muy' (t + h)emu du. (4.8)
u=O
This transformation can only be used when m > 0; otherwise we solve the problem using
an explicit method that converges to the result of transformation (4.8). Our next task is
to find m by adopting the same idea we used in finding Equation (4.7). Ifthe function
y'(t+h)emu is expanded using a Maclaurin series, we have
y'(t+h)emu = y'(t) + u[y"(t) + my'(t)] + O(h2
).
Now, we choose m such that y"(t) + my'(t) = 0 and O(s) becomes zero. From these
conditions, we have to choose m = -y"(t)/y'(t) or my'(t) = -y"(t). By considering the last
conditions and applying them to ( 4.8), we have
y(t+h) = y(t) + [y'(t)/m][l-e·mh] + [O(llm3)or O(h3
)]. (4.9)
Now, we investigate Equation (4.9) to find an appropriate explicit method used for
solving problems form::; 0. Let us take h small enough. We know perfectly well that
e-mh can be approximated using a Taylor series such as [1 - rnh + (rnhi/2 + O(h3
)]. By
substituting this result into ( 4.9), we have
y(t+h) = y(t)- y'(t)rnh2/2 + O(h3
)
= y(t) + y"(t)h2/2 + O(h\
82
This last result is a Taylor series which is equivalent to a second-order explicit RungeKutta
(RK2) method.
The numerical computation then is done as follows: If m > 0 we apply Equation
[4.9], otherwise we apply RK2. Suppose dold = f(t,y(t)) and dnew = f(t+h,y(t+h)). An
algorithm for an explicit exponential method can be given as follows:
Fori= 1 toN I* N is number of variables* I
end do
Compute doldi = f/t,y(t))
Compute dnewi = ~(t+h,y(t+h))
If doldi ::t:- 0 then
endif
Calculate mi = (doldi- dnewi)/(h*doldi)
Ifmi > 0 then
/*Calculate y(t+h) using an explicit exponential method *I
Yi(t+h) = y/t) + doldD-exp(-mi*h)]/mi
endif
if doldi = 0 or mi ~ 0 then
/*Calculate y(t+h) using RK2 */
y/t+h) = y/t) + 0.5*(doldi + dnewJ*h
endif
In order to adjust the stepsize so the computation results will compromise with the
required tolerance, we use a formula
83
(
relative error tolerance ]
112
h = *h .
new II current relative local error! I current ·
(4.10)
where llcurrent relative local errorll can be written as
N {I estimated local error for ith variable!}
m~ ' i=i IYi (t + h)l+£
and £ is a positive parameter. If IYi(t+h)l > £ then the relative error TOL is controlled,
otherwise the absolute error £TOL is controlled. The integration stepsize then is chosen
as h = min {hnew• 2hcurrent}. Usually the local error can be taken as the difference between
two methods with different orders [see Appendix D]. The local error for RK2 in this
method can be taken as the difference between the solutions obtained by RK2 and those
of the Euler method. For an exponential method, the local error is estimated as the
difference between the solutions obtained by the exponential method (4.9) and those of
Equation ( 4.11) given below
y/t+h) = Yi(t)*exp[yi'(t)h/yi(t)] + O(h2
). ( 4.11)
4.6 ODE Solvers
Several ODE solvers that are available in the public domain (see Appendix A)
will be investigated using test cases such as: Kidney problems proposed by Scott and
Watts [17, 21, 79], an autocatalytic reaction proposed by Robertson [21, 56], problem D4
of Enright [39], and a problem proposed by Gupta and Wallace [53] . Among the ODE
solvers used are MEBDF, LSODE, EPSODE, and VODE. All of the software including
the explicit methods are coded in the FORTRAN language.
84
Gear was the first person who developed a widely-used code for solving stiff
ODE problems. His early code named DIFSUB [1, 17, 44, 57, 99, 105] was developed
using BDF methods for stiff problems and Adams methods for nonstiff problems.
Among other variants of DIFSUB are STIFF [17] and GEAR (modification of STIFF
and created by Gear and Hindmarsh [17]). GEAR itself has many variants including
GEARB, GEARS, GEARBI, GEARBIL, BEARIB, GEARY, GEARST [17].
LSODE (Livermore Solver for ODEs) is a package created by Hindmarsh
[http://www.netlib.org/odepack/index.html] that is based on combinations of GEAR and
GEARB. It can be used to solve stiff or nonstiff problems. The Jacobian for the case of
stiff problems is treated as either full or banded, and as either user-supplied or internally
approximated by difference equations [ 1].
EPISODE (Experimental Package for Integration of Systems of Ordinary
Differential Equations) was created based on DIFSUB [16]. it uses implicit BDF
methods of orders one through five for stiff ODEs and Adams-Moulton methods of
orders one through twelve for nonstiff ODEs. Both BDF and Adams-Moulton methods
are in Norsieck forms. Among variants of EPISODE are EPSODE, EPISODEB, and
EPISODEIB [17]. EPSODE [http://www.netlib.org/odelindex.html] is an early version
of EPISODE modified by Byrne and Hindmarsh.
MEBDF is a package program for stiff ODEs created by Cash and Cosidine [20,
21, http://netlib.att.com/netlib/ode/index.html, http://www.netlib.org/odelindex.html]. It
is developed using modified extended BDF methods. The advantage of this approach is
85
that the methods is A-stable up to order 4 and A( a)-stable up to order 9, while pure BDF
methods themselves are A-stable only up to order 2 [28].
VODE is a package program for ODEs created by Brown, Byrne, and Hindmarsh
[ 11 , http:/ /netli b .att.com/netlib/ ode/index.html, http://www .netli b. org/ ode/index.html].
This program is a combination of EPISODE and EPISODEB. It uses variable-coefficient
Adams-Moulton and BDF methods in Nordsieck form. It treats Jacobian matrices as full
or banded.
STINT is a package program created by J.M. Tendler, T.A. Bickart, and Z. Picel
[17, http://www.netlib.org/toms/index.html]. It was the first ODE solver developed based
on cyclic composite linear multistep methods. This solver can solve ODE problems that
belong to stiffly stable.
I'
j
CHAPTER V
ANALYSIS AND DISCUSSION
The criteria of evaluation will be based on the efficiency of the codes using
parameters such as computing time (CT), number of integration steps (NS), number of
evaluations off(t,y) (NF), number of evaluations ofthe Jacobian (NJ), actual stepsize (H),
and the actual order of the methods (P). The graphical interpretations for all solutions
using methods and input test cases mentioned in section 4.1, will be compared with the
graph of solutions of the same problem solved using the NDSolve function of
Mathematica. NDSolve[ {equation_!, equation_2, ... , equation_n}, y, {t, t_min, t_max}]
is the function used in Mathematica to solve ODE problems numerically. The NDSolve
function will automatically apply the Adams predictor-corrector method if it detects that
the problem is non-stiff and will apply the backward differentiation formulae (Gear's
method) if it detects that the problem is stiff. Since NDSolve can detect the stiffness of an
ODE problem automatically, users do not need to know what stiffness is or even to be
aware of its existence. The selection method for solving stiff and non-stiff systems of
ODE for NDSolve is based on the method proposed by Petzold [95].
The modified Euler method will be implemented without the facility to calculate
the Jacobian matrix numerically; users should be able to supply the Jacobian
86
87
matrix analytically. Since the Jacobian matrix of an ODE problem is not always easy to
find, it is recommended that users use Mathematica or Maple or Matlab to get the
Jacobian aforementioned. In this thesis, the Jacobian matrix of system of ODE is
calculated using Mathematica. One example of Mathematica instructions for finding the
Jacobian matrix of problem 2 of section 4.4 is given as follows:
fly1_,y2__]=0.01-(l +(yl + 1 OOO)(y1 + 1))(0.01 +y1 +y2);
g[y1_,y2__]=0.01-(l +y2/\2)(0.01 +y1 +y2);
jacob={ {D[fly1 ,y2],yl ],D[fly1 ,y2],y2]},
{D[g[y 1 ,y2],y1 ],D[g[y1 ,y2],y2]}};
MatrixFormuacob]
Output is given as follows:
-1-(1 +y1 )(1 OOO+y1 )-(1 001 +2y1 )(0.01 +y1 +y2)
-1-y22
-1-(1 +y1 )(1 OOO+y1)
-1-y22 -2y2(0.0 1 +y 1 +y2)
5.1 Implementation ofModified Euler Method
The program for this method is given in Appendix F under the name of "Program
5". For the case of problem 1 mentioned in Chapter 4.4, the solutions are compared with
the exact solution, while for problem 2 of Chapter 4.4, the results will be compared with
the output resulted from Mathematica.
To calculate problem 1 of Chapter 4.4 from 0 to 10, this program took 520 steps,
while Lambert [77] claims that his program with the same modified Euler method only
;
88
took 501 steps and Euler's rule with the maximum allowable steplength 0.001 (the
maximum steplength that is allowable for stability of the explicit Euler method for this
problem) took 10,000 steps. The comparison between the solution of problem 1 that
resulted from this modified Euler method and the exact solution is given in Table G.5 of
Appendix G.
II)
~
>-
3-e y 1'=-y 1-0.5y2-0.5y3, y 1(0)=-1
y2'=-0.5y 1-1 000.75y2+999.25y3, y2(0)=1
y3'=-0.5y1+999.25y2-1 000.75y3, y3(0)=3
2 f\,
... y1(t)
1r~ • y2(t)
• y3(t)
0~ -~~, •••••••
-1
tAxis
Figure 5.1 Problem 1 of Section 4.4
Using Modified Euler Method
10
This "Program 5" needs about 548 steps to integrate problem 2 of section 4.4
from 0 to 10, while Lambert [77] claims that his program with the same modified Euler
89
method only needs 379 steps to solve the same cases, and a fourth-order Runge-Kutta
method took 3296 steps for a tolerance ofO.Ol (the maximum steplength that is allowable
for stability for the fourth-order Runge-Kutta method for this problem). So far, there is
no known analytic solution of this problem, so that solutions of problem 2 of section 4.4,
will be compared with the solutions from Mathematica. The comparison table resulted
from this program and Mathematica software is given in Table G.6 of Appendix G. The
graphical solutions this problem 2 produced by the modified Euler method is shown in
Figure 5.2.
0.10 I
y1'=0.01-[1+(y1 +1000)(y1 +1 )][0.01 +y1 +y2], y1(0)=0
y2'=0.01-(1+y2
2)(0.01+y1+y2). YiO)=O . ~ ...- __._. •
0.05 --i / • y1(t)
...-
• ..- • y2(t)
__. ..
1_.-
-~ 0.00
<{ >- m\..
-0.05
-0.10
? 4 6 8
tAxis '
Figure 5.2 Problem 2 of Section 4.4
Using Modified Euler Method
10
90
5.2 Implementation of an Exponential Method
This method does not need the Jacobian matrix. A method for selecting an initial
stepsize proposed by Gl~dwell et al. [ 48] has been implemented in this program. The
FORTRAN coding of this method is given in Appendix F under the name of "Program
6". The two problems used in section 5.1 are also used as examples solved with this
method.
To solve the problem 1 of section 4.4 from 0 to 10, this program took 9659 steps
for a local tolerance of 1 0-4. The comparison table between solutions of this problem 1
solved by the exponential method and the exact solution is given in Table G.7 of
Appendix G. The solutions of this problem solved with the exponential method are
shown in Figure 5.3.
~,.,
3
2
I \
y,'=-y,-0.5y2·0.5y,, y,(0)=-1
y;=-0 .5y ,-1 000. 75y2+999 .25y,, y 2(0 )= 1
y 3'=-0 .5y, +999 .2 5y2-1 000. 75y 3 , y3(0 )= 3
• y,(t)
• y 2(1) ... , r''..___ y,(t)
0
I
-1 / 1 Axis
Figure 5.3 Problem 1 of Section 4.4
Using Exponential Method
....
91
This "Program 6" takes 4369 steps to integrate problem 2 of section 4.4 from 0 to
10 with a local tolerance of 10-6
• The comparison between solutions of this problem 2
solved with the exponential method and solutions resulted from Mathematica is given in
Table 5.4. The solutions of this problem 2 solved with the exponential method is shown
in Figure 5.4.
0.10
0.05
Ul ~ 0.00
>-
-0.05
-0.10
y 1 '=0.01-[1 +(y 1+1 OO)(y 1+1 )][0.01 +y 1+y2], y 1 (0)=0
y2'=0.01-(1 +y/)(0.01 +y 1+y2), y2(0)=0 .... ----------
........ _. ..
~---..
• y1(t)
___. • Yz(t)
Figure 5.4 Problem 2 of Section 4.4
Using Exponential Method
92
5.3 Analysis of Kidney Problems
Due to Scott and Watts [20] and Byrne and Hindmarsh [16], kidney problems
with several initial conditions given in Problem 1 of section 4.1 are considered as non-stiff
and stiff problems. The graphical solutions of these problems using Mathematica
with coding instructions given in Appendix F under the name "Mathematica 1 ", are
shown in Figure 5.5 and Figure 5.6.
y y
1=0 .9902 688359 1=0.990283499
21 w.:r.t!.!.":::.~ .. :.:.: t 5
-2 0.20.40.60.8 1 t
-4 -5 0.2 0.40.60.8 1
-6
-81 -1
-101 , 1 1 F 1 F ........... -15
y y
1=0 .9925211341 1=1. 0304879856
.,,.~- , •• rtl ....
1~1 •'' _,. ..
I''' 10
, ..
,, ~~- ,. .... t'' J ... ' , ,,
t • J t
-51 0.2 0.4 0.6 0.8 1
-10
0.2 0.4 0.6 0.8 1
-15
-201 ' -20
Figure 5.5 Kidney Problems with A= 902688359, 0.0990283499,
0.9925211341, and 1.0304879856
y
1=0.99
140 tr
120 :it
100 !
80 r
60 J
40 / J
~ 20 ....... ""! ...... ' "t ......_-~-·-·.
0.2 0.4 0.6 0.8 1
y
600000
500000
400000
300000
200000
100000
1=0.
.•. !I
'I
.~
;t
f
I
.J
/ '
1~
' .... ··X: •• t
0.20.40.60.8 1
y
50000
40000
30000
20000
10000
1=0.9
.. 'I
>I
v
I
,I
I
I
i'
i
I ,J
l
~~ -+ •.. a• t
0.20.40.60.8 1
Figure 5.6 Kidney Problems with 'A= 0.99, 0.9, and 0
93
Main programs for calling MEBDF, VODE, LSODE, and EPSODE are given in
Appendix F with names "Program 1 ", "Program 2", "Program 3", and "Program 4"
respectively. All programs were run on a Sequent machine with a tolerance of 10-6
. The
performances of these four routines for the case of kidney problems are given in Table
G.l of Appendix G.
94
Byrne and Hindmarsh (17] solved this problem using LSODE for the choices of
A.= 0.99026833, 0.99, 0.9. They stated that the last two choices are stiff problems, and
the first choice is non-stiff problem. By investigating Table G.1, we can conclude that
kidney problems belong to stiff problems for the choices of A.= 0.99, 0.9, and 0.
It is difficult to analyze the performances of these four routines based on pure
CPU time of computations, because the DTIME routine supplied by Sequent can produce
different execution time for different executions of the same program with the same
problem. This occurs because Sequent uses symmetric multiprocessing (SMP) computing
concept; Multiprocessing design in which any CPU can be assigned any application task.
In all computations using MEBDF, LSODE, EPSODE, and VODE, we tried to run each
problem as many as 20 times and took the average of execution times. It is more
convenient to analyze the performances based on NS, NF, and NJ. For the case of
/....=0.9902688359, to compute this problem from 0 to 1, EPSODE needs 68 steps, 128
evaluations of f(t,y), and 23 evaluations of the Jacobian matrix of f(t,y), VODE needs 74
steps, 105 evaluations off(t,y), and 2 evaluations ofthe Jacobian matrix off(t,y), LSODE
needs 7 steps, 15 evaluations of f(t,y), and 7 evaluations of the Jacobian matrix of f(t,y),
and MEBDF needs 77 steps, 147 evaluations off(t,y), and 2 evaluations of the Jacobian
matrix of f(t,y). Since evaluations of matrices are relatively expensive, the investigation
performances will be based on the number of evaluations of the Jacobian matrix. For the
choices of A.= 0.99, 0.9, and 0.0, LSODE seems to have the best performance among the
four packages. For other choices of A., conclusions cannot be taken merely based on NJ,
since NS and NF are varying.
I
I
95
5.4 Analysis of Autocatalytic Problems
Autocatalytic reaction pathway problems are frequently called Robertson
problems to honor Robertson who was the first man who proposed these problems for
chemical kinetics problem [21, 56, 79, 1 05]. These Robertson problems will be solved
numerically using Mathematica, LSODE, MEBDF, VODE, and EPSODE, from 0 to
4. OElO. The graphical solutions will be given as functions of log10 oft. Figure 5.7
shows graphical solutions that resulted by solving these problems using the Mathematica
instructions given in Appendix F under the name "Mathematica 2".
y
1
0.8
0.6
0.4
0.2
2
r * * *-*
I
4 6
*y1
<>y2
*y3
8 10
Figure 5.7 Output ofMathematica for Robertson Problem
log10 t
96
The performances of the MEBDF, LSODE, VODE, and EPSODE routines when
solving Robertson problems are given in Table G.2 of Appendix G. Figure 5.8 is a
graphical solutions of Robertson Problems solved with EPSODE.
Figure 5.8 Robertson Problem Solved with EPSODE
The graphs that are given in Figure 5.7 and Figure 5.8 are plotted from t=IO to 4. OEIO.
Hindmarsh and Byrne [79] claimed that these problems are stiff because for the
very early stages of computation, the stepsizes should be maintained small, even though
gradually the stepsize can become larger and larger as t goes to infinity. Since NS, NF
and NJ among the four routines are varying from one to another, we cannot comment on
the best method for solving these problems among the four routines.
97
5.5 Analysis Problem D4 of Enright et al.
Enright et al. [39] proposed this problem as one example of a stiff problem. They
gave the third equation as
y3' = - 0.013 y1 - 1000 y1 Y3- 2500 Y2 Y3·
Later on, Shampine [103] commented that logically, this equation cannot match with the
initial conditions Y~> y2, and y3 at t=O that are given as 1, "1, and 0 respectively. In his
opinion, for those conditions, y3' (0) = -0.013 < 0, so that y3(t) should be negative for a
very small value of t. Shampine then improved the problem by suggesting the third
equation as
Y3' = 0.013 YI- 1000 YI Y3- 2500 Y2 Y3·
Using the Mathematica instructions given in Appendix F under the name "Mathematica
3", we get graphical solutions as shown in Figure 5.9.
y1
1
0.8
2S
y3
-6
3.S 10
+---------- t
2S so
y2
1
0.9
1
0.8
0.6
0.4
0.2
y
2S
..... t
so
""·~· ::--.---.---.. --.
_y1 -
-- y2
-Y3
-1---------- t
0 10 20 30 40 so
Figure 5.9 Problem D4 of Enright et al. Solved with Mathematica
,,;t
98
Table G.3 of Appendix G is a table of performances of the MEBDF, VODE,
LSODE, and EPSODE routines. A graphical solutions of this D4 problems solved with
MEBDF is shown in Figure 5.1 0.
1.00 --.._
. --------- ..__ ------------ ...._ ____ -...._
---...._ ---- ....__
0.75
1/)
-~
>. 0.50 • y 1(t)
• y 2(t)
0.25 ... y 3(t)
0.00
0 10 20 30 40 50
taxis
Figure 5.10 Problem D4 ofEnright et al. Solved with MEBDF
By investigating Table G.3 and assuming that evaluating the Jacobian matrix is
relatively more expensive than evaluating the function f(t,y), it can be concluded that for
these problems, MEBDF gives the best performance among the four packages.
-~·-·~- ~,J_
5.6 Analysis of Gupta and Wallace's Problem
Mathematica needs 1523 steps to solve these problems from 0 to 10. In order to
get Mathematica successfully to integrate the problems from 0 to 10, we have to supply
to the NDSolve routine of Mathematica, a parameter MaxStep with a value not less than
1523. Otherwise, Mathematica will automatically stop the execution, and gives the
message:
NDSolve :: mxst:
Maximum number of 500 steps reached at the point 3.27871.
Without supplying the parameter MaxStep, Mathematica will automatically give 500
steps as a default value for MaxStep. Figure 5.11 shows graphical solutions resulting
from Mathematica for these problems.
y1
200