- Email: [email protected]
The evaluation of numerical techniques for solution of stiff ordinary differential equations arising from chemical kinetic problems
T h e Evaluation of N u m e r i c a l Techniques for Solution of Stiff Ordinary Differential Equations Arising from Chemical Kinetic Problems D. Shyan-Shu Shieh, Y. Chang and G.R. Carmichael Chemical and Materials Engineering, Universityof Iowa, Iowa City, Iowa 52242, USA ABSTRACT Ordinary differential equations with widely scattered eigenvalues (stiff O.D.E.'s) occur often in the studies of reaction network problems. Five numerical methods, including two methods based on Backward differentiation formulas, a modified Runge-Kutta-Fehlberg method, a method based on PECE Adams formulas, and an improved semi-implicit Euler method are evaluated by comparing their performance when applied to test systems. The test systems represent different combinations of linearity and nonlinearity, small and large dimension, real and complex eigenvalues, and slightly stiff and very stiff problems. The relative merits and dificiencies of £he methods are discussed. Key words:
Stiff Ordinary Differential
Equations;
Chemical Modeling;
Introduction Ordinary differential equations are stiff when the eigenvalues of the Jacobian matrix are widely scattered. This problem was first pointed out to be :timportant in practical numerical problems by Curtiss and Hirshfelder [i] in 1952. Numerical analysts didn't focus their attention to the practical difficulty of solving stiff differential equations until ten years later when Dahlquist [2] pointed out that instability is the cause of the difficulty. Thereafter there has developed a large literature on the numerical solution of stiff systems. The need to solve stiff ordinary differential equations arises frequently in engineering and scientific analysis. Such situations occur in the study of exothermic chemical reaction in tubular reactors [3,4,5], process dynamics and control [6, 7], transient models for fluidized bed reactors [8], distillation column design [9], and reversible enzyme kinetics [I0]. Discretization of partial differential equations also gives rise to stiff differential equations. Examples include two-dimensional reactor models [ii, 12], transients for catalyst particles [13], and problems in fixed bed adsorption [14]. Other examples of stiff differential systems arising in engineering applications include circuit theory [15], nuclear reactor control rod problems [i0], impulse transmission in nerve fibres, mechanical systems, and electron distribution in lasers [16]. The need to solve stiff differential equations also arises frequently in environmental analysis. For example, nearly 90% of the computation time in simulating urban airshed and regional-scale photochemical oxidant and acid deposition problems is spent solving the stiff equations arising from the chemical processes.
Photochemistry.
[21], and Warner in 1977 [22]. However, since the time of these studies additional developments and improvements in numerical methods for O.D.E.'s have occurred. Therefore, an update to these studies was initiated and the results are reported in this paper. The objective of this work is to review some of the common techniques available for solving stiff systems of equations and to highlight the differences and similarities of these methods and the effect of these differences on their relative performances by solving various test systems with widely different characteristics. The test problems represent a wide spectrum of problems types: i.e., linearity and nonlinearity, small and large dimension, real and complex eigenvalues of Jacobian matrix, and slightly stiff to very stiff problems. The numerical results are analyzed and discussed in terms of efficiency and accuracy. This information may help users to select the most efficient method associated with the specific properties of their problems. A Review of Methods Available for Solution of Stiff Ordinary Differential Equations In this section the framework of the numerical methods available for solving stiff systems is described. Five representative methods [i.e., GEAR, EPISODE, GERK, DE AND SIEU) were selected for evaluation. The following references are recommended for a more complete description: [16] and [36] for GEAR and EP[SODE;[24] for DE.[23]for GERK,[25! for SIEU. The basic set of stiff equations represented by
to be solved is
dz d-~ = f(z' t)
(i)
The numerical solution of stiff systems of equations is difficult and time consuming (and thus, costly) because the stiffness imposes severe step-size limitations on the numerical method. Therefore, it is necessary to carefully select the technique for a given problem. However, it is difficult to identify the "best" method a priori because it will depend on the nature o ~ the problem and also on any imposed constraints (e.g., computation time, accuracy, difficulty of implementation, etc.)
where~is
There have been several survey studies conducted assessing the efficiency of the various ODE solvers and reviewing the difficulty of solving stiff problems. The first review was by Seinfeld, Lapidus, and Hwang in 1970 [17]; this was followed chronologically by Bjurel et al. in 1970 [18], Gelinas in 1972 [19], Willoughly in 1974 [i0], ~ r i g h t in 1975 [20], Schiesser and Lapidus in 1976
where Y is the solution vector, n represents the value at time t n, and k = 1 and £ = p-I for non-stiff conditions, and k = p and £ =0 for stiff conditions. In the case of EPISODE, the determination of step size h n and order p are made from step to step, as is the evaluation of ani and b . , while in GEAR h is unchanged until p+2 n3 n
Paper accepied ~ March 1988 Referee: Dr. C. Seigneur
28
Environmental Software, 1988, Vol. 3, No. 1
the vector containing
the unknowns.
EPISODE AND GEAR This implementation is based on the backward differentiation formulas k £ Y = Z Y + h Z b f (Y-Y-n-j' ) --n i=l ani -- n-I nj= 0 nj -tn-j
(2)
~) Computation~ Mechanics Pub~cations
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
c o n s e c u t i v e successful steps. In Eq. 2, bno is always greater than zero, so both routines are implicit types. Either successive substitution or N e w t o n - R a p h s o n is used to obtain the ~ values.
scheme to solve chemical reaction These problems can be written as
dt
1
t
f
n
k$ i k
P
_mI
Yi,n-1
+ h
Y. i,n =
(3)
Di
form becomes
~P£. n £ l,n-i
~Y m m,n-i
(12)
~i~Yk'n-i by
S*hn_ I
n
D.I dt
i=l ..... p
:~-
(4)
~
h
(13) n
nIi
3"=i tn-I
max{ IYi,n - Yi, n-ll/Yi, n-1 }
-tn-j
This is the A d a m s - B a s h f o r t h formula of order p, while c o r r e c t o r is the Adams-Moulton formula of order p p-i +h Z n i= 1
b f 1~-i
+hbf n o n
limit for accuracy.
The algorithm in this scheme is that Y.l,n-i. is used with hn_ I to predict Yi n and then h n is obtained from Eq. 13. With this new ~tep size h n , Y:l ,n is recalculated . by Eq. 12. This method has the feature of A-stablllty.
(5)
j#i
=Y ~-I
(10) m
The step size can be controlled
tn_ 1
t-t
i
where D 3 " is the rate coefficient of the reaction zn i. . • which Yi itself ms one of the reactants and P. Is the rate coefficient of the reaction in which Yi wlll be generated. They proposed a modified Euler method which is implicit to itself only and explicit with regard to the others; i.e., replace Yi, n-1 with Yi,n on the right hand side of Eq. Ii, thus
where s is a given
=
~
(ii)
I + h n ZDJ . i,n-I J
and
Y
l
j
+ 7.P£ ~Ym, n-i £ i,n-i
This method is based on the non-stiff formulation of Eq. 2, (i.e., k = i, £ = p-l) , and a traditional p r e d i c t o r - c o r r e c t o r procedure. Shampine and Jordon [24] s p e c i f i e d the equations of the predictor step as
=~-'-
i
which under the transformation of difference Y Y i,n - i,n-1 = ZD j ~Yk,n-i h Yi,n-i . i,n-I n ] k~i
DE
ai
problems.
dY.
The formulas in Eq. 2 belong to the family of linear m u l t i - s t e p methods which, as pointed out by D a h l q u i s t [2], only the implicit type of order 1 and 2 are Astable. Gear [26], [36], has shown that these methods are stiffly stable through order 6. Simply stated a method is defined as stable if numerical errors do not grow w i t h time. A method is A-stable if it is stable for all problems where the real parts of the eigenvalues are n e g a t i v e [2]. A method is A(u) stable if there are s t a b i l i t y restrictions on the imaginary part of the eigenvalues. Further details on stability are found in D a h l q u i s t [2] and Finnlayson [30]. Updated versions of GEAR's method are available (cf., LSODE).
P Y = Y + h Zaif-l--n --n --n-i n i=l
network
the
CHEM EQ: This technique developed by Young and Boris [32] is a second-order predictor, iterated corrector scheme which utilizes an asymptotic treatment for the stiff equations.
(6)
For equations
written
in the form
t bi = h I n [~0 Hi =
f
n tn-1
H i dt
(7)
t - tn_ j tn_i -tn_ j
dY. ___~i= p.(y) dt i --
(8)
_ D.(Y)Y. I -- 1
The stiff equations
Yi,n = The values of a i and b i are computed as the s t e p size and order are adjusted. Error estimation in this package is b a s e d on the imbedding technique as E
=
n
Y
-
n,p+l
Y
to predictor:
(15)
2 Ti,n_ 1 + ~t
and corrector:
m
(9)
n,p
are solved according
Yi,n_l[21"i,n_l - 6t] + 26t Xi,n_ 1 P i , n - 1
•
j%i
(14)
Y. l,n
=
{
6t m-i ~ - [Ti
m-1 +
ri,n-l]
[Pi
+
Pi,n-i ]
(16)
where Y is obtained by the procedure u s i n g a n,p+1 predictor of order p and a corrector of order p+l, Y nwp • . . is o b t a l n e d by the procedure of uslng both p r e d l c t o r and corrector of order p. This method is c o n d i t i o n a l l y stable.
+ Yi,n-l[~7 -1 + Ti,n-i 6t]}/
[T m-1 i
+ r i , n + 6t]
with r m. ~ I/D~, and m is the iteration number. The time step t~ advance from the n-i time level to the nth level can be selected according to the relation
GERK This m e t h o d is a modified explicit fourth-order "~ Runge-Kutta formula developed by Fehlberg [23] in order to minimize the truncation error. W a t t and Shampine [24] c o m p l e m e n t e d it by applying the Richardson-type technique to control the s t e p size. D a h l q u i s t [2] has shown that no e x p l i c i t method is A-stable and thus, G E R K is a conditionally stable method. SIEU Originally,
Preussner
and Brand
[25] presented
this
6 t = eimin
Yi,n-i { Pi,n-1 - Di,n-lYi,n-i
}
(17)
~SA: This technique developed by Hesstvedt, Hov and Isaksen [33] is based on expressing the problem in the form of Eq. 14. If the production (i.e., Pi ) and the destruction terms (D i) are assumed constant over the time
Environmental Software, 1988, Vol. 3, No. 1
29
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
step 6t, i.e.,
then this e q u a t i o n can be solved analytically,
Pi Yn = ~ ( l_eDi 6t) + Yn-I e - D i 6 t 1
become complex. The r e a c t i o n m e c h a n i s m and the initial c o n c e n t r a t i o n s are given in Tables 1 and 2, respectively. The e i g e n v a l u e s at two d i f f e r e n t times, t=0 and t = 1 . 0 x l 0 -5, are listed in Table 3.
(18) Table
The accuracy of this method is dependent on the careful choice of St. Hesstvedt et al. [33] p r e s e n t a scheme by w h i c h Eq. 18 is used only for species where ~ t/10 < T. < I00 A t. For species where T. >> ~ t a -1 simple e x p l i c i t l E ~ l e r ' s i n t e g r a t i o n is recommended, and for s p e c i e s where r. ~ t/10 values are c a l c u l a t e d based on the P a s s u m p t i o n of local equilibrium, i.e., y =__nn n D n The last three methods discussed, i.e., SIEU, C H E M EQ, and QSSA all r e p r e s e n t techniques which make use of asymptotic or "quasi-linear" assumptions. D e s c r i p t i o n of the Test Systems The n u m e r i c a l techniques d i s c u s s e d in the last s e c t i o n are a p p l i e d to three test p r o b l e m s in order to d e m o n s t a t e their applicability. These e x a m p l e s cover broad ranges of the following characteristics: (I) s y s t e m dimension, (2) stiffness ratio, (3) real or c o m p l e x e i g e n v a l u e s of the Jacobian matrix, and (4) l i n e a r i t y or nonlinearity. For these tests only the S I E U m e t h o d from the g r o u p of a s y m p t o p i c methods is d i s c u s s e d .
I
C h e m i c a l Reactions and Rate C o n s t a n t s for P r o b l e m II
chemical kinetics rate c o n s t a n t ......................................................... I. 2. 3. 4. 5. 6. 7. 8. 9. i0. 11. 12. 13. 14. 15. 16. 17.
03 + h~ ÷ 0('D) + 02 0 ('D) + H20 ~ 20H 0 ('D) + M + O + M O + 02 + M + 03 + M S02 + OH ÷ HS03 SO2 + HO2 + OH + SO 3 H02 + NO + OH + NO2 ON + NO ÷ HNO2 HNO2 + hv + OH + NO NO2 + hv ÷ NO + 0 CO + OH ÷ CO2 + H HO2 + H02 ÷ H202 + 02 2HNO2 + NO + NO2 +H20 NO + NO2 + H20 ÷ 2HN02 H + 02 + M ÷ H02 + M H02 + OH + H20 + 02 NO + 03 + NO2 + 02
KI=3.4E-3 K2=7.35E 5 K3=I.20E 5 K4-1.70E-5 K5=I00E 3 K6=1.20 K7=3.00E 2 K8=3.00E 3 K9=3.00E-2 KI0=0.35 K11=2.00E 2 KI2-8.40E 3 K13=0°I0 KI4=2.00E-6 KI5-I.70E-3 KI6=3.00E 5 KI7=2.00E I
P R O B L E M I: N O N L I N E A R S Y S T E M OF R E A C T I O N R A T E EQUATIONSI'? [17] Table dYl/dt=-kl +k2*Y2*Y3 d Y 2 / d t = k I *YI - k 2 * Y 2 * Y 3 - k 3 * Y 2 * * 2 dY3/dt=k3*Y2**2 I.C.:
(19) xl
(0)
=
Y2 Y3
(0) = 0 (0) = 0
k2*Y3
k2*Y2
A =kl
-k2*Y3-2*k3*Y2
-k2*Y2
0
2*k3*Y2
0
Initial C o n c e n t r a t i o n s ppm for Problem II
i. 2. 3. 4. 5. 6. 7. 8. 9.
5 320 1.00E-If 1.00E-2 3.70E-5 1.40E 4 6.5E-2 1.00E-2 4.61
i0. 11. 12. 13. 14. 15. 16. 17.
Table
[-B+(B**2-4*C)**0.5]/2
i = 2,3
w h e r e B = 0.04 + 1.0x104*y3 + 0 . 6 x 1 0 8 * y 2 C =
0.24xlO7*y2 + 0.60xlO12*y2**2
w h e r e at t = 0,
~
= 0
X
I as t ÷ ~
,
Y1 = 0, ~ = 0 I
= 0,
X
2 Y2 X
2
~ -0.04 3
= O, = 0,
with kl = 0.04, k3=3.0x107, k2=variable. ratio (SR) is i n c r e a s e d as k2 increases. becomes very stiff as time increases.
Y3 A
3
= 1 -4 - -ix I 0
The s t i f f n e s s Also the s y s t e m
3.32 1.00E-11 6.60E-16 2°00E 5 1.49 2.90E-7 3.20 1.00E-4
3
The E i g e n v a l u e s of the M a t r i x A for Problem II
I. 2. 3. 4. ~. 6. 7. 8. 9. 10. 11. 12. 13.
t=0 -0.I06E12 -0.272E 9 -0.272E 7 -0.181E 5 -0.131E
4
0.422E 3 -0.895E 2 -0.460E 2 -0.468E 1 -0.174E 0 -0.303E-3 -0.327E-5 -0.162E-6
t=l. 0E-5 -0.I06E12 -0.272E 9 - 0.272E 7 -0.181E 5 0.132E
SMOG
This system is a s i m p l i f i e d c h e m i s t r y model s t u d i e d by C a r m i c h a e l and Peters [27]. It c o n s i s t s of 17 reactions and 17 species, and r e p r e s e n t s a nonlinear, m o d e r a t e l y large-dimension, very stiff problem. Initially, the eigenvalues of the Jacobian m a t r i x are real, and r i g h t after i n t e g r a t i o n begins, some o f t h e m
Environmental Software, 1988, Vol. 3, No. 1
4
-0.131E 4 -0.101E 3 O + (0.125E i) -0.681E 0 - (0.125E i) 0.117E-I - 0 . 4 0 2 E - 4 + (0.585E-4) - 0 . 4 0 2 E - 4 - (0.585E-4) -0.185E-8 -0.681E
All the r e m a i n i n g e i g e n v a l u e s P R O B L E M II: R E A C T I O N N E T W O R K FOR P H O T O C H E M I C A L
30
NO 2 0 0('D) 02 03 OH SO 2 SO 3
= 0
~. = 1
and and
CO CO 2 B BNO 2 HO 2 H 20 H202 HSO 3 NO
(20)
which %s singular, with ki=0.04, k2=l.0xl04, k3 1.0xl0-. ~he three e i g e n v a l u e s of the m a t r i x A are I
in
i
This system represents one that is nonlinear, with small d i m e n s i o n and real e i g e n v a l u e s of the J a c o b i a n matrix. Here, YI, Y2, Y3 are the mole fractions of the r e a c t i n g species. The J a c o b i a n m a t r i x A of the system is given by -kl
2
are zero
i
i i i
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al LINEAR SYSTEM W I T H E A S Y - T O - A D J U S T P R O B L E M III: E I G E N V A L U E S [28] d Y 1 / d t = -10*YI + a ' Y 2 dY2/dt= -~'YI - 10*Y2 d Y 3 / d t = -4"Y3 (21) d Y 4 / d t = -Y4 dY5/dt = -0.5"Y5 d Y 6 / d t = -0. I*Y6 w i t h I.C. : Yi(0)=l i=i,2,...,6 and Eigenvalues: - 0.I, -0.5, -i, -4, and -I0+
(a) ¢ -i for i = I, 2, ...,6
The s o l u t i o n
to Eq. 21 is
Y l = e x p ( - 1 0 t ) * s i n ( u t) (22)
Y2=exp(-10t)*cos(ut) and
Yi=exp
i= 3, 4, 5, 6
(-~it) ,
time c o n s u m i n g m a t r i x operations such as LU decomposition, multiplication, or inversion, with TOL = 10 -4, d i v i s i o n by zero occured in DE. That is the reason why T O L = I 0 - ~ was set. It should also be noted that b e c a u s e the roundoff c o n s t a n t of double p r e c i s i o n (which contributes to part of the p r o p a g a t i o n error) is smaller, c a l c u l a t i o n with double p r e c i s i o n always runs faster than with single p r e c i s i o n . Also, it is found that c a l c u l a t i o n with single p r e c i s i o n often restricts the a p p l i c a b i l i t y of the methods. This later phenomenon has been p o i n t e d out by Chan [29] in 1977. In order to see the effect of stiffness ratio on the p e r f o r m a n c e of the methods, the p a r a m e t e r k2 in Problem I is reset to 108 instead of i0^4. Under these c o n d i t i o n s the s t i f f n e s s ratio is increased by a factor of l04 and the results are given in Tables 6 and 7. As shown, G E R K and DE were unable to provide a solution within the cutoff time period (180 CPU-s) and that EPISODE and SIEU p r o v i d e d the more accurate solutions. However, as the tolerance is d e c r e a s e d to 10 -8 the GEAR model supplies very a c c u r a t e results at about a tenth of the CPU of EPISODE for this problem.
where 8 3 = -4 84 -I
85 86
-0.5 =
-0. I
The m a j o r reason in c h o o s i n g i m a g i n a r y part of c o m p l e x ~ can be value and it is e x p e c t e d that when e i g e n v a l u e s close to the imaginary methods will be handicapped (i.e., backward differentiation).
this system is that the a d j u s t e d to be any the J a c o b i a n has >~ axis c e r t a i n n u m e r i c a l those based on
TEST PROBLEM R E S U L T S S o l u t i o n s to the three test problems using the n u m e r i c a l methods d e s c r i b e d were o b t a i n e d using a PRIME 750 system and e x c e p t when noted all testing was done in d o u b l e precision. F u r t h e r m o r e it must be noted that both G E A R and E P I S O D E have stiff and n o n s t i f f routines and four d i f f e r e n t ways for e s t i m a t i n g the J a c o b i a n matrix. In this c o m p a r i s o n work the J a c o b i a n m a t r i x of each test s y s t e m was e s t i m a t e d analytically, and the number "21" b e h i n d the code name r e p r e s e n t s the stiff routine w h i l e "ii" r e p r e s e n t s the nonstiff routine. A n u m b e r of factors might be c o n s i d e r e d in a c c e s s i n g the r e l a t i v e merits of the d i f f e r e n t methods. The c o m p a r i s o n s p r e s e n t e d here are r e s t r i c t e d to the factor of cost (CPU) and accuracy. However, the reader may d e c i d e to w e i g h t the d i f f e r e n t factors in a p a r t i c u l a r case a c c o r d i n g to their p a r t i c u l a r a p p l i c a t i o n and needs. The f o l l o w i n g statistics are compiled for each test: CPU TIME: F L N EVAL: J A C O B EVAL: LAST H: LAST ORDER: STEPS TAKEN: TOL:
Total c o m p u t a t i o n time required to solve the p r o b l e m The number of f u n c t i o n e v a l u a t i o n s required to solve the p r o b l e m Total number of J a c o b i a n m a t r i x e v a l u a t i o n s The last step size used at the e n d p o i n t of the i n t e g r a t i o n The last order of method used at the e n d p o i n t of the i n t e g r a t i o n The number of steps r e q u i r e d to solve the problem S p e c i f i e d error tolerance. In SIEU, we s p e c i f i e d s (See Eq. 13) instead of TOL.
T a b l e s 4 and 5 r e p r e s e n t the c o m p u t a t i o n results for the n o n l i n e a r problem I, which is of d i m e n s i o n three, and w i t h real eigenvalues. %~e e x a c t solutions at t=10 are reported by M i c h e l s e n [28] as Yi-0.841370, Y 2 = 0 . 1 6 2 3 4 " i 0 -4, Y3=0.158614. It can be seen from t h e s e d e t a i l e d results that GEAR21 takes the least c o m p u t i n g time w h i l e GERK has the most a c c u r a t e solution. Due to the s t a b i l i t y limit, GERK and DE have the s m a l l e s t last s t e p size and take the most c o m p u t i n g time although both codes do not need the Jacobian m a t r i x e v a l u a t i o n and the
P r o b l e m II r e p r e s e n t s the s t i f f e s t case in our test systems. ~ stiffness ratio (SR) in this problem is as large as i0 . F i n l a y s o n has p o i n t e d 3 o u t that "typi~ally, SR=20 is not stiff, SR=I0 is stiff, and SR=I0 is very stiff" [30]. Moreover, this system represents a nonlinear, modestly large-dimension (consisting of 17 d e p e n d e n t variables) problem. The ~calculated results for this problem are summarized in Tables 8 and 9. The integrations for GEAR11 with T O L > 0 . 0 1 and E P I S O D E 1 1 with TOL < 1.0 (which are the n o n s t i f f routines) and DE and GERK were halted as they took more than 180 seconds CPU time. This system is a p r o b l e m that does not have an analytic solution for comparison. However, it is known that a severe error tolerance would lead to accurate solution. Meanwhile, the c o n c e n t r a t i o n s obtained by d i f f e r e n t codes are about the same order of magnitude. Strictly speaking, GEAR and EPISODE gave a l m o s t the same solution and those of SIEU had up to 40% d e v i a t i o n in some species. Again, it is noted that there is neither error e s t i m a t i o n nor t o l e r a n c e c o n t r o l in this i m p l e m e n t a t i o n of the SIEU method. C o n s i d e r i n g the factor of cost only in Problem II, SIEU p e r f o r m e d best. For methods based on BDF (background d i f f e r e n c e s formula), a matrix inversion is required. Usually, the computing time of the matrix i n v e r s i o n is p r o p o r t i o n a l to the second to third power of the d i m e n s i o n of the test system. Problem III represents a linear ODE set with e i g e n v a l u e s of the J a c o b i a n matrix close to the imaginary axis. The p e r f o r m a n c e of the methods are summarized in Tables i0 and ii. Due to the roundoff constant of the computer the final solution of Y1 and Y2 always d e v i a t e d from the exact solutions and the solutions of Y3 to Y6 are used to j u s t i f y if the results are right or not. As seen, G E A R and SIEU failed to give the right solutions. On the other hand, EPISODE, GERK, and DE all offered the same solutions, and cost is the only factor needed to be considered. For this s l i g h t l y stiff p r o b l e m (SR=I02) DE and GERK w o r k e d better (because they do not require the time c o n s u m i n g m a t r i x operation) than EPISODE. In addition, the d e t e r i o r a t i o n of the p e r f o r m a n c e of EPISODE21 is caused by the poor s t a b i l i t y regions of the higher order formulas it used. For example, we ran this problem with = 1.0, i0, i00 r e s p e c t i v e l y (results are p r e s e n t e d e l s e w h e r e [31]) and it was found that in all cases, DE was the most e f f i c i e n t and GERK and EPISODE11 were of the same order of efficiency. It is also interesting to note that the s t a b i l i t y limit is no longer of concern in this p r o b l e m so that the four s u c c e s s f u l routines have the same order of m a g n i t u d e of last step size.
Environmental Software, 1988, Vol. 3, No. 1
31
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al |
Table Results
Numerical
TOL
CPU TIME
LAST H
4 for
Problem
LAST ORDER
I
STEPS TAKEN
FCN EVAL
JACOB EVAL
(SEC) .......................................................................................................................... GEAR ii GEAR 21 EPISODE ii EPISODE 21 GERK DE SIEU
1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-6 1.0E-3 2.0E-3 3.0E-3
0.42 0.41 0.43 0.39 1.96 1.03 0.98 0.80 33.1 46.9 5.42 1.91 0.92
1.3 2.7 I.I 2.3 4.1E-2 0.39 1.2 2.4 1.5E-3 7.5E-4 2.0E-2 4.1E-2 7.5E-2
3 2 3 1 2 2 4 3 4 5 1 1 1
46 40 51 39 373 119 143 114 6714 14591 7501 2352 913
62 47 67 44 440 150 194 139 114860 31181 15002 4704 1826
14 13 15 13 II0 87 97 91 -
Table 5 P e r c e n t Errors on S o l u t i o n Yi for Problem I at t=10
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
V
.
TOL
GEAR
1.0E-4
ii
1.0E-2
GEAR
1.0E-4
21
1.0E-2
EPISODE
.
1.0E-2
EPISODE
1.0E-4
21
1.0E-2
GERK DE
1.0E-4 1.0E-6
SIEU
1.0E-3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2.0E-3 3.0E-3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Y2*I04
0.84136 0.0010% 0.86418 2.7710% 0.84149 0.0152% 0.84983 1.0064% 0.84136 0.0002% 0.84127 0.010% 0.84135 0.002% 0.84050 0.103% 0.84137 0.84111 0.030% 0.84169 0.0385% 0.84844 0.8406% 0.87789 4.3416%
1.0E-4
ii
.
Y1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Y3
0.15689 3.3578% 0.63011 288.459% 0.16244 0.0622% 0.17008 4.7678% 0.16234 0.000% 0.16185 0.302% 0.16234 0.000% 0.16095 0.856% 0.16234 0.15088 7.059% 0.16227 0.04131% 0.16127 0.6591% 0.15455 4.7986%
0.15862 0.0057% 0.13575 14.4104% 0.15848 0.0807% 0.15014 5.3394% 0.15861 0.001% 0.15871 0.063% 0.15862 0.009% 0.15948 0.547% 0.15861 0.15886 0.160% 0.15883 0.1362% 0.16229 2.3176% 0.18083 14.0063%
Table 6 N u m e r i c a l Results for P r o b l e m I with K 2 = I 0 8 .......................................................................................................................... TOL
CPU LAST LAST STEPS FCN JACOB TIME H ORDER TAKEN EVAL EVAL (SEC) .......................................................................................................................... GEAR ii GEAR 21 EPISODE 11 EPISODE 21 GERK DE SIEU
32
10 -4
i.ii
10 -4
0.82
10 -4 10 -4 10 -4 10 -4 0.001 0.0005
2
157
474
86
0.13
1
113
226
62
9.34
0.0054
2
2716
2943
248
1.47
0.81
4
229
369
Iii
0.139 0.0157
1 1
7134 14078
14268 28156
-
5.68 10.62
"~ 2.9
Environmental Software, 1988, Vol. 3, No. 1
The .Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
Table Solution with
Yi
for
K2=I08
7 Problem at
I
t=10
...................................................................................................... TOL
Y1
Y2*I06
Y3*I03
....................................................................................................... GEAR ii
10 - 4
1.0076 (0.128%)
1.95581 (156%)
-0.7584 (245%)
GEAR 21
10 - 4
1.0014 (0.0066%
-2.76930 (463%)
-0.1358 (126%)
EPISODE ii
10 -4
0.99948 0%
0.76326 ~ 0%
0.52357 ~ 0%
EPISODE 21
10 -4
0.99948 0%
0.76323 (0.004%)
0.52359 (0.004%)
GERK
10 -4
-
- -
DE
10 -4
-
- -
SIEU
0.001
0.99982 (0.034%)
0.0005
0.99967 (0.019%)
~
0.4572 (40%)
0.8748 (6.71%)
0.7506 (1.65%)
0.5327 (1.74%)
Table 8 N u m e r i c a l R e s u l t s for P r o b l e m II at t = 5 0 0
TOL
GEAR II GEAR 21 EPISODE ii EPISODE 21 SIEU
CPU TIME (SEC)
i. 0 E - 2 1.0 I. 0 E - 4 io 0 E - 2 1,0 1.0
5.03 3.27 9.42 5.12 3.81 5.84
I. 0 E - 4 i. 0 E - 2 1.0 0.01 0.05 0.i0 1.00
12.22 8.33 5.82 15 • 63 3.92 2.42 1.12
LAST H
LAST ORDER
STEPS TAKEN
3.1E 1 1.2E 2 4.5E 1 9.5E 1 4.0E 2 1.9E 2
2 2 4 3 1 2
116 54 285 124 65 75
5.0E 1.4E 2.2E i. 44 7.75 1.62E 5.5E
4 3 2 1 1 1 1
301 158 77 3487 709 362 51
1 2 2
1 2
FCN EVAL
JACOB EVAL
150 60 327 152 80 85
28 19 49 30 25 60
449 222 89 6974 1418 724 102
97 87 68 -
Environmental Software, 1988, Vol. 3, No. 1
33
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al Table 9
Numerical Solution Yi for Problem II at t=500
TOL
GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
Y1
Y2*10 -2
Y3*I0 II
Y4
Y5*106
Y6*10 -4
0.01 1.00 10 -4 0.01 1.00 1.00
2.6702 2.6660 2.6702 2.6708 2.7401 2.6715
3.2233 3.2233 3.2233 3.2233 3.2226 3.2233
1.7375 1.7212 1.7232 1.7242 1.6551 1.7236
1.4884 1.4784 1.4884 1.4884 1.4198 1.4907
9.4037 9.6524 9.4182 9.4089 10.1278 9.4255
1.3998 1.3998 1.3998 1.3998 1.3998 1.3998
10 -4 0.01 1.00 0.01 0.05 0.I0 1.00
2.6702 2.6719 2.6719 2.7732 2.7684 2.7720 2.7037
3.2233 3.2233 3.2233 3.2224 3.2224 3.2227 3.2255
1.7232 1.7230 1.7242 1.3946 1.4232 1.4012 1.7538
1.4884 1.4880 1.4938 1.6219 1.6213 1.6227 1.5996
9.4186 9.4373 9.4532 5.4623 5.7192 5.5736 9.1437
1.3998 1.3998 1.3998 1.3998 1.3998 1.3998 1.3996
TOL
Y7*I0 -2
Y8
Y9
YI0
Yll,107
Y12,1015
..................................................................................................... GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
0.01 1.00 10 -4 0.01 1.00 1.00
6.5209 6.5100 6.5209 6.5210 6.5392 6.5215
3o06~6 3.0067 ~ 3.0673 3.0577 3.0081 3.0863
1.6492 1.6191 1.6490 1.6447 1.7419 1.6593
4.8024 4.8424 4.8026 4.8069 4.7783 4.7900
6.1793 6.2317 6.1804 6.1860 6.1492 6.1643
1,6300 1,6754 1,6301 1,6233 1o8636 1.6361
10 -4 0.01 1.00 0.01 0.05 0.i0 1.00
6.5209 6.5211 6.5227 6.5171 6.5164 6.5182 6.5374
3.0672 3.0584 3.1241 3.1338 3.0588 3.2313 4.2285
1.6490 1.6468 1.6783 2.3254 2.2505 2.3030 1.8368
4.8026 4.8052 4.7679 4.0503 4.1349 4.0967 5.0598
6.1805 6.1838 6.1358 5.2098 5.3197 5.2691 6.5336
1.6302 1,6372 1.6509 0.9726 1.0267 0,9943 1.5454
YI4*I02
YI5*I06
YI6*I0
YI7*I0
TOL
GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
34
Y13"I0-5
0.01 1.00 10 -4 0,01 1.00 1.00
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
5.0960 5.2375 5.0958 5.0748 5.8259 5.1146
8.7794 8.7990 8.7767 8.7815 8.6570 8.7923
1.3800 2.0011 1.3834 1.4799 1.9797 1.1932
4.4427 3.3321 4.439L 4.4043 4.0353 4.4675
10 -4 0.01 1.00 0.01 0.05 0.i0 1.00
1.9999 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
5.0964 5.1182 5.1610 3.0436 3.2123 3.1087 4.8445
8.7769 8.7755 8.8447 6.8001 6.9881 6.8661 8.1412
1.3840 1.4732 0.8152 1.7607 1.6760 1.8483 2.4357
4.4374 4.3936 4.5437 4.5987 4.3991 4.8281 7.7700
Environmental Software, 1988, Vol. 3, No. 1
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al COMPARISON OF QSSA AND GEARS METHOD FOR A REALISTIC AIR POLLUTION APPLICATION The analysis of air pollution chemistry problems, such as photochemical oxidant and acid rain formation, requires the solution of large sets of very stiff equations. For example, the analysis of 50-100 chemical species involved in one-to-several hundred reactions is not uncommon for these problems. Studies of chemical mechanisms typically make use of one of the Gear's methods(i.e., Eq. 2). However, when multidimensional analysis of the problems is performed, the computational burden of using Gear's method becomes too great, and one of the asymptotic methods is used. To better understand the trade-offs between the asymptotic methods and Gear's method a multiday box model calculation using the Lurman, Loyd, Atkinson [34] chemical mechanism (referred to hereafter in this paper as LLA) was performed. The LLA mechanism consists of 53 chemical species and 112 chemical reactions. For comparison purposes the Q6SA method as implemented in the STEM-II model (cf, Carmichael, et al., [35]) is compared with results obtained using the Gear's method. The QSSA method is as indicated by Eq. 18. However, the time steps are fixed at 30 seconds during the daytime , 60 seconds at night time, and 12 seconds at sunrise and sunset. The comparisons between Gear's method and the QSSA method are performed for the gas phase kinetics in a box model. In order to emphasize only the chemistry, the :~ physical processes such as emissions, depositions, dilutions, transports, etc. are not included. A test simulation adopted initial concentrations typical of polluted areas near sunrise: 40 ppb of NO, i0 ppb of NO 2, 20 ppb of S02, 1 ppb of 03, 30 ppb of alkanes, 8 ppb of ethene, 7 ppb of alkenes, I0 ppb of aromatics, 7 ppb of aldehydes and 1 ppb of ketones. Typical background concentrations of CH 4 and CO of 1650 and 125 ppb, respectively, were also included throughout this simulation. The ambient temperature and relative humidity were held to be fixed at 20 C and 50%, respectively. A multi-day simulation for the summer solstice at 40 N latitude was employed. The comparisons for key species are shown in Figure i. The percentage relative error based on results from the Gear's method is listed in Table 12. Except for relatively large discrepancies in species like PAN, the QSSA method predicts reliable results when compared with Gear's method. Figure 2 shows the percent relative error in the mass balance of sulfur and nitrogen. The use of the QSSA method often fails to maintain the mass balance. Moreover, the strong coupling between more reactive species result in numerical instability if a relatively long time step is used. Implementation of a linear transformation (or lumping) technique allows the circumvention of these two main difficulties. These procedures are explained in detail in Hesstvedt et ai.[33]. Results using the QSSA method with lumping are presented in Table 12 and Figures 1 ~ 2. Clearly, the lumping technique maintains the mass exactly for sulfur and within 0.5% for nitrogen throughout the simulation, except at the peak occuring around 18:00 of the first day. This seems to arise from numerical stiffness associated with a dramatic change at sunset. On the other hand, for the case without lumping, the strong discrepancies are seen in nitrogen compounds, which are..relatively more reactive than sulfur compounds. The QSSA method executes ~ l0 times faster than Gear's method for this application.
suited for specific types of applications. In order to assess and quantify the attributes of these techniques, five representative methods (i.e., GEAR, EPISODE, GERK, DE, and SIEU) were applied to test problems. These test problems were chosen to span the range of problem types, including: linear and nonlinear; small and large dimension; real and complex eignevalues; and slightly stiff and very stiff conditions. Results of the numerical test calculations indicate that: i.
If the accuracy of solutions is considered alone, GERK was the best candidate for nonstiff problems. Due to the small stability limit, this is not recommended to be used for stiff ordinary differential equations. DE was quite suitable for use in solving moderately stiff systems (SR < 102) with eigenvalues of Jacobian matrix close to the imaginary axis. For very stiff systems, SIEU was the best method at less stringent tolerance because of its short computation time. It worked better than GEAR and EPISODE if the dimension of problem is greater than ~ 15 - 20. However, it did not work as well as GEAR and EPISODE at more stringent tolerances. The two BDF methods, GEAR and EPISODE were roughly comparable in terms of accuracy and computing time, and each resulted in significant savings of time over the explicit methods DE and GERK. EPISODE resulted in the higher accuracy in each example and GEAR worked more efficiently, except in obtaining solutions under conditions when the eigenvalues were close to the imaginary axis. If it is necessary to choose one as an all-purpose stiff solver, EPISODE with routine 21 is recommended. For use in large air pollution models the asymptotic methods offer a good blend of accuracy and conservation properties of the solution.
2.
3.
4.
5. 6.
ACKNOWLEDGEMENT Partial support of this work was provided by California Air Resources Board through a subcontract with Sigma Research. REFERENCES [i]
Curtiss, C.F. and Hirschfelder, J.O. Integration of Stiff Equations, Proc. of National Adad. Sci., U.S.A., 1952, 38, 3, 235.
[2]
Dahlquist, G. G., "A Special Stability Problem for Linear Multistep Methods", BIT 3, pp. 27-43 (1963).
[3]
Eschenroeder, A. Q., Boyer, D. W., and Hall, J.G., "Nonequilibrium Expansions of Air with Coupled Chemical Reactions", Phys. Fluids, 5, NO. 5, pp. 615-624 (1962).
[4]
Amundson, N. R., "Some Further Observations on Tubular Reactor Stability", Can. J. Chem. Eng., 43, 2, 49-55 (1965).
[5]
[6]
Davison, E. J., "The Numerical Solution of Large Systems of Linear Differential Equations", AICh. E. J., 14, l, 46 (1968).
[7]
Mah, R.S.H., Michaelson, S. and Sargent, R.W.H., "Dynamic Behavior of Multi-component and Multistage Systems", Chem. ~ . Sci., I__~7,619-639 (1962).
[8]
Luss, D. and Amundson, N.R., "Stability of Batch Catalytic Fluidized Beds", A I C h . E . J . , 14, 211 (1968).
[9]
Slavickova, A., Hlavecek, V., and Kubicek, M., "Symposium on Computer in the Design and Erection of Chemical Plants", Karlory Vary (1975).
SUMMARY Many techniques are available for use in the numerical solution of stiff ordinary differential equations. However, each of the different techniques has its own relative strengths and weaknesses and is best
Amundson, N. R., Luss, D., "Qualitative and Quantitative Observations on the Tubular Reactor", Can. J. Chem. Eng., 46, 6, 425-433 (1968).
Environmental Software, 1988, Vol. 3, No. 1
35
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al [i0]
Willoughby, R.A., Ed., "Stiff Differential Systems", the IBM Research Symposia Series, Plenum Press, New York, NY (1974).
[ii]
[29]
Chan, Y. N. I., et al., "Solution of Stiff Differential Equations and the Use of Imbedding Techniques", Ind. En~. Chem. Fundam., 17, 3, 133 (1978).
Finlayson, B. A., "Packed Bed Reactor Analysis by Orthogonal Collocation", Chem. En~. Sci., 26, 7, 1081 (1971).
[30]
Finlayson, B.A., "Nonlinear Analysis in Chemical Engineering", McGraw-Hill, Inc. (1980).
(12]
Froment, G. F., "Fixed Bed Catalytic Reactors", Ind. En~. Chem., 59___, 2, 18 (1967).
[31]
[13]
Lee, J.C.M. and Luss, D., "The Effect of Lewis Numer on the Stability of a Catalytic Reaction", AICh. E. J., 16, 620 (1970).
Shieh, D. S-S, "The Evaluation of Numerical Techniques for Solution of Stiff Ordinary Differential Equations Arising from Chemical Kinetic Problems", M.S. Thesis, Chemical and Materials Engineering, University of Iowa (1983).
[32] [14]
Fleck, R.D., Kirwan, D.J. and Hall, K.R., "Mixed Resistance Diffusion Kinetics in Fixed-Bed Adsorption under Constant Pettern Conditions", Ind. En~. Chem. Fundam., 12, i, 95 (1973).
Young, T.R. and Boris, J.P., A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reaction flow problems, J. Phys. Chem., 1977, 81, 2424.
[33] [15]
Brayton, R.K., Gustavson, F.G. and Liniger, W., "A Numerical Analysis of the Transient Behavior of a Transistor Circuit", IBM J. Res. Develop., I0___, 4, 292 (1966).
[16]
Byrne, G.D. et al., "A Comparison of Two ODE codes: GEAR and EPISODE", Com~ and Chem. E., 1,133-147 (1977).
Hesstvedt, E., Hov, O., and Isaksen, I.S.A., Qusisteady-state approximations in air pollution modeling: comparison of two numerical schemes for oxidant prediction, Int. J. of Chem. Kinetics, 1978, x, 971. Lurmann, F.W., Lloyd, A.C. and Atkinson, R., A chemical mechanism for use in long-range transport/acid deposition computer modeling, JGR, 1986, 91, 10905.
[17]
Beinfeld, J.H., Lapidus, L. and Hwang, M., "Review of Numerical Integration Techniques for Stiff Ordinary Differential Equations", Ind. Eng. Chem, Fundam., 9, 2, 266 (1970).
[34]
[18]
Bjurel, G., Dahlquist, G.G., Lindberg, B., and Oden, S., "Survey of Stiff Ordinary Differential Equations", Report No. NA70-11, Department of Information Processing, The Royal Inst. of Tech., Stockholm, Sweden (1970).
[19]
Gelinas, R.J., "Stiff Systems of Kinetic Equations", J. of Comput. Physics, 9__, 222-236 (1972).
[20]
Enright, W. H., Hull, T. E. and Lindberg, B., "Comparing Numerical Methods for Stiff Systems of ODEs", Bit 15, pp. 10-48 (1975).
[21]
Lapidus, L. and Schiesser, W.E., Ed., "Numerical Methods for Differential Systems", Academic Press, Inc. New York (1976).
[22]
Warner, D. D., "The Numerical Solution of the Equations of Chemical Kinetics", J. Phys. Chem., 81, 25, 2329 (1977).
[23]
Fehlberg, E., "Low-order Classical Runge-Kutta Formulas with Stepsize Control and Their Application to Some Heat Transfer Problems" NASA TR R-315, July (1969).
[24]
Shampine, L. F. and Gordon, M. K., "Computer Solution of Ordinary Differential Equations", W. H. Freeman and Company (1975).
[25]
Preussner, P. R. and Brand, K. P., "Application of a Semi-implicit Euler Method to Mass Action Kinetics", Chem. En~. Sci., 36___,I0, 1633 (1981).
[26]
Gear, C. W., "Numerical Initial Value Problems in Ordinary Differential Equations", Prentice-Hall,F Englewood Cliffs, NJ, pp. 158-166 (1971).
[27]
Carmichael, G. R. and Peters, L. K., "A Regional Model for SO x Transport in the Troposphere", Atmospheric Environment, 18___,937 (1984).
[28]
Michelsen, M. L., "An Efficient General Purpose Method for the Integration of Stiff Ordinary Differential Equations", A I C h . E . J . , 22, 3, 594 {1976).
36
Environmental Software, 1988, Vol. 3, No. 1
[35]
Carmichael, G.R., Peters, L.K. and Kitada, T., A second generation model for regional-scale transport~chemistry~deposition, Atmos. Environ., 1986, 20, 173.
[36]
Gear, C.W., The automatic integration of ordinary differential equations, Cumm. Acm, 1971, 14, 185.
.~
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh eI al
=• .,.o
LEgENO
''iiiiiiiiiii
,,
L
°i
°0 _.e~
O * - MSIM 14/ LUMPING o - M S I M W/O LUMPINg ~.- gERR'S METHOD
./r
~f=
o 12
DRY
18
2t
6
12
18
24
ORT 2 TIME O" CLOCK
1
OR~( 3
i
12
IB
OR~
24
6
12
1B
24
ORY 2 I|ME O'CLOCK
1
ORY 3
y................ :I m @' SO= z
c~
c-I LEGEND + - MSIM W/ LUMPING o - MSIM W/O LUMPINg - gERR'S M~THOD
z
8
,
~
~ ,.=,@]
~ ~
LEGEND + - MSIM W/ LUMPINg o - MSIM W/O LUMPING
"
_,.4_=,~ ~,~,"; ;- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
¢0 ~_ .3 E. 6
,
,
,
,
,
,
,
12
18
24
6
12
18
24
DRY 1
OR'~ 2 T ] ME O' CLOCK
....
S
,2
DRY 3
,~
DRY
=,
1 TIME
Figure
I.
Comparison between Gear's method.
the Q S S A method
~
~
~
~'~
ORY 2 O' CLOCK
DRY
and
N ! TROGEN
z
el: _J ,0--
tr~
o o . H l TH LUHP ING
_
z
u~
w~
o
~
xx H I T H O U T
SULFUR
LUMPING
............. .,~
td >
......
J
.........
? e
,'2
,;
DRY l
Figure
2.
2'4
;
,'~
,'8
fi
DAY 2 T I'M[':O' CLOCK
Comparison of the mass balance with and without lumping.
s
DR1 3
-
Environmental Software, 1988, Vol. 3, No. 1
37
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
Table I0 Numerical Results for Problem III w i t h u = 1.0 ....................................................................................................... TOL
GEAR ii GEAR 21 EPISODE 11 EPISODE GERK DE SIED
CPU TIME
LAST H
LAST ORDER
STEPS TAKEN
FCN EVAL
JACOB EVAL
1.0E-6
1.99
0.75
3
202
269
25
1.0E-6
1.57
0.83
5
159
191
24
1.0E-6
11.67
0.026
8
943
1510
144
1.0E-6 1.0E-6 1.0E-6 0.5 0.6
14.49 10.37 7.30 2.37 1.25
0.012 0.021 0.021 0.092 0.122
5 4(5) I0 1 1
1901 2787 927 907 289
2524 20798 1874 1814 578
227
Table II N u m e r i c a l S o l u t i o n Yi for Problem III at t=20 w i t h u = 1 ...................................................................................................... Y1
Y2
Y3
Y4*I09
Y5*I05
Y6*I0
...................................................................................................... GEAR Ii GEAR 21 EPISODE ii EPISODE 21 GERK DE SIEU (s=0.5) SIEU (s=0.6)
2.26E-8
1.46E-8
-2.311E-20
1.69
4.55
1.35
-1.13E-f1
-3.48E-II
2.16E-II
5.82
4.59
1.35
1.83E-87
-6.99E-88
1.805E-35
2.06
4.54
1.35
1.83E-87
-6.98E-88
1.805E-35
2.06
4.54
1.35
1.83E-87 1.83E-87 3.96E-66
-6.99E-88 -6.99E-88 -1.34E-66
1.805E-35 1.805E-35 2.683E-31
2.06 2.06 4.10
4.54 4.54 5.41
1.35 1.35 1.36
2.29E-60
6.82E-61
7.031E-30
5.40
5.83
1.37
1.26E-87
5.65E-88
1.850E-35
2.06
4.54
1.35
SOL'N
Table 12 C o m p a r i s o n b e t w e e n the QSSA Method and Gear's M e t h o d
Species
% Relative After 1 day MSIM(1) MSIM(2)
HNO 3 SO 2 Sulfate 03 ALKA ETHE AROM HCHO H202 PAN KET
0.24 -0.03 -0.62 -0.91 -0.24 0.49 -0.55 -8.22 -9.71 -11.5 -1.80
*
**
38
R e l a t i v e Error With Lumping
=
Gear
-
-5.96 0.ii -1.83 -2.50 -0.39 2.60 2.05 -9.45 -10.2 "~-15.6 -3.80
MSDIM)/Gear*I00
Environmental Software, 1988, Vol. 3, No. 1
Error* After 2 days MSIM(1)** (MSIM(2)*** 0.46 -0.16 0.45 -0.81 -0.59 0.39 -2.96 -6.80 -2.52 -13.3 -1.31
-5.70 0.i0 -0.75 2.57 0.53 4.79 1.73 -7.97 -3.67 -17.8 -3.30
Stiff Ordinary Differential
Equations;
Chemical Modeling;
Introduction Ordinary differential equations are stiff when the eigenvalues of the Jacobian matrix are widely scattered. This problem was first pointed out to be :timportant in practical numerical problems by Curtiss and Hirshfelder [i] in 1952. Numerical analysts didn't focus their attention to the practical difficulty of solving stiff differential equations until ten years later when Dahlquist [2] pointed out that instability is the cause of the difficulty. Thereafter there has developed a large literature on the numerical solution of stiff systems. The need to solve stiff ordinary differential equations arises frequently in engineering and scientific analysis. Such situations occur in the study of exothermic chemical reaction in tubular reactors [3,4,5], process dynamics and control [6, 7], transient models for fluidized bed reactors [8], distillation column design [9], and reversible enzyme kinetics [I0]. Discretization of partial differential equations also gives rise to stiff differential equations. Examples include two-dimensional reactor models [ii, 12], transients for catalyst particles [13], and problems in fixed bed adsorption [14]. Other examples of stiff differential systems arising in engineering applications include circuit theory [15], nuclear reactor control rod problems [i0], impulse transmission in nerve fibres, mechanical systems, and electron distribution in lasers [16]. The need to solve stiff differential equations also arises frequently in environmental analysis. For example, nearly 90% of the computation time in simulating urban airshed and regional-scale photochemical oxidant and acid deposition problems is spent solving the stiff equations arising from the chemical processes.
Photochemistry.
[21], and Warner in 1977 [22]. However, since the time of these studies additional developments and improvements in numerical methods for O.D.E.'s have occurred. Therefore, an update to these studies was initiated and the results are reported in this paper. The objective of this work is to review some of the common techniques available for solving stiff systems of equations and to highlight the differences and similarities of these methods and the effect of these differences on their relative performances by solving various test systems with widely different characteristics. The test problems represent a wide spectrum of problems types: i.e., linearity and nonlinearity, small and large dimension, real and complex eigenvalues of Jacobian matrix, and slightly stiff to very stiff problems. The numerical results are analyzed and discussed in terms of efficiency and accuracy. This information may help users to select the most efficient method associated with the specific properties of their problems. A Review of Methods Available for Solution of Stiff Ordinary Differential Equations In this section the framework of the numerical methods available for solving stiff systems is described. Five representative methods [i.e., GEAR, EPISODE, GERK, DE AND SIEU) were selected for evaluation. The following references are recommended for a more complete description: [16] and [36] for GEAR and EP[SODE;[24] for DE.[23]for GERK,[25! for SIEU. The basic set of stiff equations represented by
to be solved is
dz d-~ = f(z' t)
(i)
The numerical solution of stiff systems of equations is difficult and time consuming (and thus, costly) because the stiffness imposes severe step-size limitations on the numerical method. Therefore, it is necessary to carefully select the technique for a given problem. However, it is difficult to identify the "best" method a priori because it will depend on the nature o ~ the problem and also on any imposed constraints (e.g., computation time, accuracy, difficulty of implementation, etc.)
where~is
There have been several survey studies conducted assessing the efficiency of the various ODE solvers and reviewing the difficulty of solving stiff problems. The first review was by Seinfeld, Lapidus, and Hwang in 1970 [17]; this was followed chronologically by Bjurel et al. in 1970 [18], Gelinas in 1972 [19], Willoughly in 1974 [i0], ~ r i g h t in 1975 [20], Schiesser and Lapidus in 1976
where Y is the solution vector, n represents the value at time t n, and k = 1 and £ = p-I for non-stiff conditions, and k = p and £ =0 for stiff conditions. In the case of EPISODE, the determination of step size h n and order p are made from step to step, as is the evaluation of ani and b . , while in GEAR h is unchanged until p+2 n3 n
Paper accepied ~ March 1988 Referee: Dr. C. Seigneur
28
Environmental Software, 1988, Vol. 3, No. 1
the vector containing
the unknowns.
EPISODE AND GEAR This implementation is based on the backward differentiation formulas k £ Y = Z Y + h Z b f (Y-Y-n-j' ) --n i=l ani -- n-I nj= 0 nj -tn-j
(2)
~) Computation~ Mechanics Pub~cations
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
c o n s e c u t i v e successful steps. In Eq. 2, bno is always greater than zero, so both routines are implicit types. Either successive substitution or N e w t o n - R a p h s o n is used to obtain the ~ values.
scheme to solve chemical reaction These problems can be written as
dt
1
t
f
n
k$ i k
P
_mI
Yi,n-1
+ h
Y. i,n =
(3)
Di
form becomes
~P£. n £ l,n-i
~Y m m,n-i
(12)
~i~Yk'n-i by
S*hn_ I
n
D.I dt
i=l ..... p
:~-
(4)
~
h
(13) n
nIi
3"=i tn-I
max{ IYi,n - Yi, n-ll/Yi, n-1 }
-tn-j
This is the A d a m s - B a s h f o r t h formula of order p, while c o r r e c t o r is the Adams-Moulton formula of order p p-i +h Z n i= 1
b f 1~-i
+hbf n o n
limit for accuracy.
The algorithm in this scheme is that Y.l,n-i. is used with hn_ I to predict Yi n and then h n is obtained from Eq. 13. With this new ~tep size h n , Y:l ,n is recalculated . by Eq. 12. This method has the feature of A-stablllty.
(5)
j#i
=Y ~-I
(10) m
The step size can be controlled
tn_ 1
t-t
i
where D 3 " is the rate coefficient of the reaction zn i. . • which Yi itself ms one of the reactants and P. Is the rate coefficient of the reaction in which Yi wlll be generated. They proposed a modified Euler method which is implicit to itself only and explicit with regard to the others; i.e., replace Yi, n-1 with Yi,n on the right hand side of Eq. Ii, thus
where s is a given
=
~
(ii)
I + h n ZDJ . i,n-I J
and
Y
l
j
+ 7.P£ ~Ym, n-i £ i,n-i
This method is based on the non-stiff formulation of Eq. 2, (i.e., k = i, £ = p-l) , and a traditional p r e d i c t o r - c o r r e c t o r procedure. Shampine and Jordon [24] s p e c i f i e d the equations of the predictor step as
=~-'-
i
which under the transformation of difference Y Y i,n - i,n-1 = ZD j ~Yk,n-i h Yi,n-i . i,n-I n ] k~i
DE
ai
problems.
dY.
The formulas in Eq. 2 belong to the family of linear m u l t i - s t e p methods which, as pointed out by D a h l q u i s t [2], only the implicit type of order 1 and 2 are Astable. Gear [26], [36], has shown that these methods are stiffly stable through order 6. Simply stated a method is defined as stable if numerical errors do not grow w i t h time. A method is A-stable if it is stable for all problems where the real parts of the eigenvalues are n e g a t i v e [2]. A method is A(u) stable if there are s t a b i l i t y restrictions on the imaginary part of the eigenvalues. Further details on stability are found in D a h l q u i s t [2] and Finnlayson [30]. Updated versions of GEAR's method are available (cf., LSODE).
P Y = Y + h Zaif-l--n --n --n-i n i=l
network
the
CHEM EQ: This technique developed by Young and Boris [32] is a second-order predictor, iterated corrector scheme which utilizes an asymptotic treatment for the stiff equations.
(6)
For equations
written
in the form
t bi = h I n [~0 Hi =
f
n tn-1
H i dt
(7)
t - tn_ j tn_i -tn_ j
dY. ___~i= p.(y) dt i --
(8)
_ D.(Y)Y. I -- 1
The stiff equations
Yi,n = The values of a i and b i are computed as the s t e p size and order are adjusted. Error estimation in this package is b a s e d on the imbedding technique as E
=
n
Y
-
n,p+l
Y
to predictor:
(15)
2 Ti,n_ 1 + ~t
and corrector:
m
(9)
n,p
are solved according
Yi,n_l[21"i,n_l - 6t] + 26t Xi,n_ 1 P i , n - 1
•
j%i
(14)
Y. l,n
=
{
6t m-i ~ - [Ti
m-1 +
ri,n-l]
[Pi
+
Pi,n-i ]
(16)
where Y is obtained by the procedure u s i n g a n,p+1 predictor of order p and a corrector of order p+l, Y nwp • . . is o b t a l n e d by the procedure of uslng both p r e d l c t o r and corrector of order p. This method is c o n d i t i o n a l l y stable.
+ Yi,n-l[~7 -1 + Ti,n-i 6t]}/
[T m-1 i
+ r i , n + 6t]
with r m. ~ I/D~, and m is the iteration number. The time step t~ advance from the n-i time level to the nth level can be selected according to the relation
GERK This m e t h o d is a modified explicit fourth-order "~ Runge-Kutta formula developed by Fehlberg [23] in order to minimize the truncation error. W a t t and Shampine [24] c o m p l e m e n t e d it by applying the Richardson-type technique to control the s t e p size. D a h l q u i s t [2] has shown that no e x p l i c i t method is A-stable and thus, G E R K is a conditionally stable method. SIEU Originally,
Preussner
and Brand
[25] presented
this
6 t = eimin
Yi,n-i { Pi,n-1 - Di,n-lYi,n-i
}
(17)
~SA: This technique developed by Hesstvedt, Hov and Isaksen [33] is based on expressing the problem in the form of Eq. 14. If the production (i.e., Pi ) and the destruction terms (D i) are assumed constant over the time
Environmental Software, 1988, Vol. 3, No. 1
29
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
step 6t, i.e.,
then this e q u a t i o n can be solved analytically,
Pi Yn = ~ ( l_eDi 6t) + Yn-I e - D i 6 t 1
become complex. The r e a c t i o n m e c h a n i s m and the initial c o n c e n t r a t i o n s are given in Tables 1 and 2, respectively. The e i g e n v a l u e s at two d i f f e r e n t times, t=0 and t = 1 . 0 x l 0 -5, are listed in Table 3.
(18) Table
The accuracy of this method is dependent on the careful choice of St. Hesstvedt et al. [33] p r e s e n t a scheme by w h i c h Eq. 18 is used only for species where ~ t/10 < T. < I00 A t. For species where T. >> ~ t a -1 simple e x p l i c i t l E ~ l e r ' s i n t e g r a t i o n is recommended, and for s p e c i e s where r. ~ t/10 values are c a l c u l a t e d based on the P a s s u m p t i o n of local equilibrium, i.e., y =__nn n D n The last three methods discussed, i.e., SIEU, C H E M EQ, and QSSA all r e p r e s e n t techniques which make use of asymptotic or "quasi-linear" assumptions. D e s c r i p t i o n of the Test Systems The n u m e r i c a l techniques d i s c u s s e d in the last s e c t i o n are a p p l i e d to three test p r o b l e m s in order to d e m o n s t a t e their applicability. These e x a m p l e s cover broad ranges of the following characteristics: (I) s y s t e m dimension, (2) stiffness ratio, (3) real or c o m p l e x e i g e n v a l u e s of the Jacobian matrix, and (4) l i n e a r i t y or nonlinearity. For these tests only the S I E U m e t h o d from the g r o u p of a s y m p t o p i c methods is d i s c u s s e d .
I
C h e m i c a l Reactions and Rate C o n s t a n t s for P r o b l e m II
chemical kinetics rate c o n s t a n t ......................................................... I. 2. 3. 4. 5. 6. 7. 8. 9. i0. 11. 12. 13. 14. 15. 16. 17.
03 + h~ ÷ 0('D) + 02 0 ('D) + H20 ~ 20H 0 ('D) + M + O + M O + 02 + M + 03 + M S02 + OH ÷ HS03 SO2 + HO2 + OH + SO 3 H02 + NO + OH + NO2 ON + NO ÷ HNO2 HNO2 + hv + OH + NO NO2 + hv ÷ NO + 0 CO + OH ÷ CO2 + H HO2 + H02 ÷ H202 + 02 2HNO2 + NO + NO2 +H20 NO + NO2 + H20 ÷ 2HN02 H + 02 + M ÷ H02 + M H02 + OH + H20 + 02 NO + 03 + NO2 + 02
KI=3.4E-3 K2=7.35E 5 K3=I.20E 5 K4-1.70E-5 K5=I00E 3 K6=1.20 K7=3.00E 2 K8=3.00E 3 K9=3.00E-2 KI0=0.35 K11=2.00E 2 KI2-8.40E 3 K13=0°I0 KI4=2.00E-6 KI5-I.70E-3 KI6=3.00E 5 KI7=2.00E I
P R O B L E M I: N O N L I N E A R S Y S T E M OF R E A C T I O N R A T E EQUATIONSI'? [17] Table dYl/dt=-kl +k2*Y2*Y3 d Y 2 / d t = k I *YI - k 2 * Y 2 * Y 3 - k 3 * Y 2 * * 2 dY3/dt=k3*Y2**2 I.C.:
(19) xl
(0)
=
Y2 Y3
(0) = 0 (0) = 0
k2*Y3
k2*Y2
A =kl
-k2*Y3-2*k3*Y2
-k2*Y2
0
2*k3*Y2
0
Initial C o n c e n t r a t i o n s ppm for Problem II
i. 2. 3. 4. 5. 6. 7. 8. 9.
5 320 1.00E-If 1.00E-2 3.70E-5 1.40E 4 6.5E-2 1.00E-2 4.61
i0. 11. 12. 13. 14. 15. 16. 17.
Table
[-B+(B**2-4*C)**0.5]/2
i = 2,3
w h e r e B = 0.04 + 1.0x104*y3 + 0 . 6 x 1 0 8 * y 2 C =
0.24xlO7*y2 + 0.60xlO12*y2**2
w h e r e at t = 0,
~
= 0
X
I as t ÷ ~
,
Y1 = 0, ~ = 0 I
= 0,
X
2 Y2 X
2
~ -0.04 3
= O, = 0,
with kl = 0.04, k3=3.0x107, k2=variable. ratio (SR) is i n c r e a s e d as k2 increases. becomes very stiff as time increases.
Y3 A
3
= 1 -4 - -ix I 0
The s t i f f n e s s Also the s y s t e m
3.32 1.00E-11 6.60E-16 2°00E 5 1.49 2.90E-7 3.20 1.00E-4
3
The E i g e n v a l u e s of the M a t r i x A for Problem II
I. 2. 3. 4. ~. 6. 7. 8. 9. 10. 11. 12. 13.
t=0 -0.I06E12 -0.272E 9 -0.272E 7 -0.181E 5 -0.131E
4
0.422E 3 -0.895E 2 -0.460E 2 -0.468E 1 -0.174E 0 -0.303E-3 -0.327E-5 -0.162E-6
t=l. 0E-5 -0.I06E12 -0.272E 9 - 0.272E 7 -0.181E 5 0.132E
SMOG
This system is a s i m p l i f i e d c h e m i s t r y model s t u d i e d by C a r m i c h a e l and Peters [27]. It c o n s i s t s of 17 reactions and 17 species, and r e p r e s e n t s a nonlinear, m o d e r a t e l y large-dimension, very stiff problem. Initially, the eigenvalues of the Jacobian m a t r i x are real, and r i g h t after i n t e g r a t i o n begins, some o f t h e m
Environmental Software, 1988, Vol. 3, No. 1
4
-0.131E 4 -0.101E 3 O + (0.125E i) -0.681E 0 - (0.125E i) 0.117E-I - 0 . 4 0 2 E - 4 + (0.585E-4) - 0 . 4 0 2 E - 4 - (0.585E-4) -0.185E-8 -0.681E
All the r e m a i n i n g e i g e n v a l u e s P R O B L E M II: R E A C T I O N N E T W O R K FOR P H O T O C H E M I C A L
30
NO 2 0 0('D) 02 03 OH SO 2 SO 3
= 0
~. = 1
and and
CO CO 2 B BNO 2 HO 2 H 20 H202 HSO 3 NO
(20)
which %s singular, with ki=0.04, k2=l.0xl04, k3 1.0xl0-. ~he three e i g e n v a l u e s of the m a t r i x A are I
in
i
This system represents one that is nonlinear, with small d i m e n s i o n and real e i g e n v a l u e s of the J a c o b i a n matrix. Here, YI, Y2, Y3 are the mole fractions of the r e a c t i n g species. The J a c o b i a n m a t r i x A of the system is given by -kl
2
are zero
i
i i i
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al LINEAR SYSTEM W I T H E A S Y - T O - A D J U S T P R O B L E M III: E I G E N V A L U E S [28] d Y 1 / d t = -10*YI + a ' Y 2 dY2/dt= -~'YI - 10*Y2 d Y 3 / d t = -4"Y3 (21) d Y 4 / d t = -Y4 dY5/dt = -0.5"Y5 d Y 6 / d t = -0. I*Y6 w i t h I.C. : Yi(0)=l i=i,2,...,6 and Eigenvalues: - 0.I, -0.5, -i, -4, and -I0+
(a) ¢ -i for i = I, 2, ...,6
The s o l u t i o n
to Eq. 21 is
Y l = e x p ( - 1 0 t ) * s i n ( u t) (22)
Y2=exp(-10t)*cos(ut) and
Yi=exp
i= 3, 4, 5, 6
(-~it) ,
time c o n s u m i n g m a t r i x operations such as LU decomposition, multiplication, or inversion, with TOL = 10 -4, d i v i s i o n by zero occured in DE. That is the reason why T O L = I 0 - ~ was set. It should also be noted that b e c a u s e the roundoff c o n s t a n t of double p r e c i s i o n (which contributes to part of the p r o p a g a t i o n error) is smaller, c a l c u l a t i o n with double p r e c i s i o n always runs faster than with single p r e c i s i o n . Also, it is found that c a l c u l a t i o n with single p r e c i s i o n often restricts the a p p l i c a b i l i t y of the methods. This later phenomenon has been p o i n t e d out by Chan [29] in 1977. In order to see the effect of stiffness ratio on the p e r f o r m a n c e of the methods, the p a r a m e t e r k2 in Problem I is reset to 108 instead of i0^4. Under these c o n d i t i o n s the s t i f f n e s s ratio is increased by a factor of l04 and the results are given in Tables 6 and 7. As shown, G E R K and DE were unable to provide a solution within the cutoff time period (180 CPU-s) and that EPISODE and SIEU p r o v i d e d the more accurate solutions. However, as the tolerance is d e c r e a s e d to 10 -8 the GEAR model supplies very a c c u r a t e results at about a tenth of the CPU of EPISODE for this problem.
where 8 3 = -4 84 -I
85 86
-0.5 =
-0. I
The m a j o r reason in c h o o s i n g i m a g i n a r y part of c o m p l e x ~ can be value and it is e x p e c t e d that when e i g e n v a l u e s close to the imaginary methods will be handicapped (i.e., backward differentiation).
this system is that the a d j u s t e d to be any the J a c o b i a n has >~ axis c e r t a i n n u m e r i c a l those based on
TEST PROBLEM R E S U L T S S o l u t i o n s to the three test problems using the n u m e r i c a l methods d e s c r i b e d were o b t a i n e d using a PRIME 750 system and e x c e p t when noted all testing was done in d o u b l e precision. F u r t h e r m o r e it must be noted that both G E A R and E P I S O D E have stiff and n o n s t i f f routines and four d i f f e r e n t ways for e s t i m a t i n g the J a c o b i a n matrix. In this c o m p a r i s o n work the J a c o b i a n m a t r i x of each test s y s t e m was e s t i m a t e d analytically, and the number "21" b e h i n d the code name r e p r e s e n t s the stiff routine w h i l e "ii" r e p r e s e n t s the nonstiff routine. A n u m b e r of factors might be c o n s i d e r e d in a c c e s s i n g the r e l a t i v e merits of the d i f f e r e n t methods. The c o m p a r i s o n s p r e s e n t e d here are r e s t r i c t e d to the factor of cost (CPU) and accuracy. However, the reader may d e c i d e to w e i g h t the d i f f e r e n t factors in a p a r t i c u l a r case a c c o r d i n g to their p a r t i c u l a r a p p l i c a t i o n and needs. The f o l l o w i n g statistics are compiled for each test: CPU TIME: F L N EVAL: J A C O B EVAL: LAST H: LAST ORDER: STEPS TAKEN: TOL:
Total c o m p u t a t i o n time required to solve the p r o b l e m The number of f u n c t i o n e v a l u a t i o n s required to solve the p r o b l e m Total number of J a c o b i a n m a t r i x e v a l u a t i o n s The last step size used at the e n d p o i n t of the i n t e g r a t i o n The last order of method used at the e n d p o i n t of the i n t e g r a t i o n The number of steps r e q u i r e d to solve the problem S p e c i f i e d error tolerance. In SIEU, we s p e c i f i e d s (See Eq. 13) instead of TOL.
T a b l e s 4 and 5 r e p r e s e n t the c o m p u t a t i o n results for the n o n l i n e a r problem I, which is of d i m e n s i o n three, and w i t h real eigenvalues. %~e e x a c t solutions at t=10 are reported by M i c h e l s e n [28] as Yi-0.841370, Y 2 = 0 . 1 6 2 3 4 " i 0 -4, Y3=0.158614. It can be seen from t h e s e d e t a i l e d results that GEAR21 takes the least c o m p u t i n g time w h i l e GERK has the most a c c u r a t e solution. Due to the s t a b i l i t y limit, GERK and DE have the s m a l l e s t last s t e p size and take the most c o m p u t i n g time although both codes do not need the Jacobian m a t r i x e v a l u a t i o n and the
P r o b l e m II r e p r e s e n t s the s t i f f e s t case in our test systems. ~ stiffness ratio (SR) in this problem is as large as i0 . F i n l a y s o n has p o i n t e d 3 o u t that "typi~ally, SR=20 is not stiff, SR=I0 is stiff, and SR=I0 is very stiff" [30]. Moreover, this system represents a nonlinear, modestly large-dimension (consisting of 17 d e p e n d e n t variables) problem. The ~calculated results for this problem are summarized in Tables 8 and 9. The integrations for GEAR11 with T O L > 0 . 0 1 and E P I S O D E 1 1 with TOL < 1.0 (which are the n o n s t i f f routines) and DE and GERK were halted as they took more than 180 seconds CPU time. This system is a p r o b l e m that does not have an analytic solution for comparison. However, it is known that a severe error tolerance would lead to accurate solution. Meanwhile, the c o n c e n t r a t i o n s obtained by d i f f e r e n t codes are about the same order of magnitude. Strictly speaking, GEAR and EPISODE gave a l m o s t the same solution and those of SIEU had up to 40% d e v i a t i o n in some species. Again, it is noted that there is neither error e s t i m a t i o n nor t o l e r a n c e c o n t r o l in this i m p l e m e n t a t i o n of the SIEU method. C o n s i d e r i n g the factor of cost only in Problem II, SIEU p e r f o r m e d best. For methods based on BDF (background d i f f e r e n c e s formula), a matrix inversion is required. Usually, the computing time of the matrix i n v e r s i o n is p r o p o r t i o n a l to the second to third power of the d i m e n s i o n of the test system. Problem III represents a linear ODE set with e i g e n v a l u e s of the J a c o b i a n matrix close to the imaginary axis. The p e r f o r m a n c e of the methods are summarized in Tables i0 and ii. Due to the roundoff constant of the computer the final solution of Y1 and Y2 always d e v i a t e d from the exact solutions and the solutions of Y3 to Y6 are used to j u s t i f y if the results are right or not. As seen, G E A R and SIEU failed to give the right solutions. On the other hand, EPISODE, GERK, and DE all offered the same solutions, and cost is the only factor needed to be considered. For this s l i g h t l y stiff p r o b l e m (SR=I02) DE and GERK w o r k e d better (because they do not require the time c o n s u m i n g m a t r i x operation) than EPISODE. In addition, the d e t e r i o r a t i o n of the p e r f o r m a n c e of EPISODE21 is caused by the poor s t a b i l i t y regions of the higher order formulas it used. For example, we ran this problem with = 1.0, i0, i00 r e s p e c t i v e l y (results are p r e s e n t e d e l s e w h e r e [31]) and it was found that in all cases, DE was the most e f f i c i e n t and GERK and EPISODE11 were of the same order of efficiency. It is also interesting to note that the s t a b i l i t y limit is no longer of concern in this p r o b l e m so that the four s u c c e s s f u l routines have the same order of m a g n i t u d e of last step size.
Environmental Software, 1988, Vol. 3, No. 1
31
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al |
Table Results
Numerical
TOL
CPU TIME
LAST H
4 for
Problem
LAST ORDER
I
STEPS TAKEN
FCN EVAL
JACOB EVAL
(SEC) .......................................................................................................................... GEAR ii GEAR 21 EPISODE ii EPISODE 21 GERK DE SIEU
1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-2 1.0E-4 1.0E-6 1.0E-3 2.0E-3 3.0E-3
0.42 0.41 0.43 0.39 1.96 1.03 0.98 0.80 33.1 46.9 5.42 1.91 0.92
1.3 2.7 I.I 2.3 4.1E-2 0.39 1.2 2.4 1.5E-3 7.5E-4 2.0E-2 4.1E-2 7.5E-2
3 2 3 1 2 2 4 3 4 5 1 1 1
46 40 51 39 373 119 143 114 6714 14591 7501 2352 913
62 47 67 44 440 150 194 139 114860 31181 15002 4704 1826
14 13 15 13 II0 87 97 91 -
Table 5 P e r c e n t Errors on S o l u t i o n Yi for Problem I at t=10
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
V
.
TOL
GEAR
1.0E-4
ii
1.0E-2
GEAR
1.0E-4
21
1.0E-2
EPISODE
.
1.0E-2
EPISODE
1.0E-4
21
1.0E-2
GERK DE
1.0E-4 1.0E-6
SIEU
1.0E-3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2.0E-3 3.0E-3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Y2*I04
0.84136 0.0010% 0.86418 2.7710% 0.84149 0.0152% 0.84983 1.0064% 0.84136 0.0002% 0.84127 0.010% 0.84135 0.002% 0.84050 0.103% 0.84137 0.84111 0.030% 0.84169 0.0385% 0.84844 0.8406% 0.87789 4.3416%
1.0E-4
ii
.
Y1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Y3
0.15689 3.3578% 0.63011 288.459% 0.16244 0.0622% 0.17008 4.7678% 0.16234 0.000% 0.16185 0.302% 0.16234 0.000% 0.16095 0.856% 0.16234 0.15088 7.059% 0.16227 0.04131% 0.16127 0.6591% 0.15455 4.7986%
0.15862 0.0057% 0.13575 14.4104% 0.15848 0.0807% 0.15014 5.3394% 0.15861 0.001% 0.15871 0.063% 0.15862 0.009% 0.15948 0.547% 0.15861 0.15886 0.160% 0.15883 0.1362% 0.16229 2.3176% 0.18083 14.0063%
Table 6 N u m e r i c a l Results for P r o b l e m I with K 2 = I 0 8 .......................................................................................................................... TOL
CPU LAST LAST STEPS FCN JACOB TIME H ORDER TAKEN EVAL EVAL (SEC) .......................................................................................................................... GEAR ii GEAR 21 EPISODE 11 EPISODE 21 GERK DE SIEU
32
10 -4
i.ii
10 -4
0.82
10 -4 10 -4 10 -4 10 -4 0.001 0.0005
2
157
474
86
0.13
1
113
226
62
9.34
0.0054
2
2716
2943
248
1.47
0.81
4
229
369
Iii
0.139 0.0157
1 1
7134 14078
14268 28156
-
5.68 10.62
"~ 2.9
Environmental Software, 1988, Vol. 3, No. 1
The .Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
Table Solution with
Yi
for
K2=I08
7 Problem at
I
t=10
...................................................................................................... TOL
Y1
Y2*I06
Y3*I03
....................................................................................................... GEAR ii
10 - 4
1.0076 (0.128%)
1.95581 (156%)
-0.7584 (245%)
GEAR 21
10 - 4
1.0014 (0.0066%
-2.76930 (463%)
-0.1358 (126%)
EPISODE ii
10 -4
0.99948 0%
0.76326 ~ 0%
0.52357 ~ 0%
EPISODE 21
10 -4
0.99948 0%
0.76323 (0.004%)
0.52359 (0.004%)
GERK
10 -4
-
- -
DE
10 -4
-
- -
SIEU
0.001
0.99982 (0.034%)
0.0005
0.99967 (0.019%)
~
0.4572 (40%)
0.8748 (6.71%)
0.7506 (1.65%)
0.5327 (1.74%)
Table 8 N u m e r i c a l R e s u l t s for P r o b l e m II at t = 5 0 0
TOL
GEAR II GEAR 21 EPISODE ii EPISODE 21 SIEU
CPU TIME (SEC)
i. 0 E - 2 1.0 I. 0 E - 4 io 0 E - 2 1,0 1.0
5.03 3.27 9.42 5.12 3.81 5.84
I. 0 E - 4 i. 0 E - 2 1.0 0.01 0.05 0.i0 1.00
12.22 8.33 5.82 15 • 63 3.92 2.42 1.12
LAST H
LAST ORDER
STEPS TAKEN
3.1E 1 1.2E 2 4.5E 1 9.5E 1 4.0E 2 1.9E 2
2 2 4 3 1 2
116 54 285 124 65 75
5.0E 1.4E 2.2E i. 44 7.75 1.62E 5.5E
4 3 2 1 1 1 1
301 158 77 3487 709 362 51
1 2 2
1 2
FCN EVAL
JACOB EVAL
150 60 327 152 80 85
28 19 49 30 25 60
449 222 89 6974 1418 724 102
97 87 68 -
Environmental Software, 1988, Vol. 3, No. 1
33
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al Table 9
Numerical Solution Yi for Problem II at t=500
TOL
GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
Y1
Y2*10 -2
Y3*I0 II
Y4
Y5*106
Y6*10 -4
0.01 1.00 10 -4 0.01 1.00 1.00
2.6702 2.6660 2.6702 2.6708 2.7401 2.6715
3.2233 3.2233 3.2233 3.2233 3.2226 3.2233
1.7375 1.7212 1.7232 1.7242 1.6551 1.7236
1.4884 1.4784 1.4884 1.4884 1.4198 1.4907
9.4037 9.6524 9.4182 9.4089 10.1278 9.4255
1.3998 1.3998 1.3998 1.3998 1.3998 1.3998
10 -4 0.01 1.00 0.01 0.05 0.I0 1.00
2.6702 2.6719 2.6719 2.7732 2.7684 2.7720 2.7037
3.2233 3.2233 3.2233 3.2224 3.2224 3.2227 3.2255
1.7232 1.7230 1.7242 1.3946 1.4232 1.4012 1.7538
1.4884 1.4880 1.4938 1.6219 1.6213 1.6227 1.5996
9.4186 9.4373 9.4532 5.4623 5.7192 5.5736 9.1437
1.3998 1.3998 1.3998 1.3998 1.3998 1.3998 1.3996
TOL
Y7*I0 -2
Y8
Y9
YI0
Yll,107
Y12,1015
..................................................................................................... GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
0.01 1.00 10 -4 0.01 1.00 1.00
6.5209 6.5100 6.5209 6.5210 6.5392 6.5215
3o06~6 3.0067 ~ 3.0673 3.0577 3.0081 3.0863
1.6492 1.6191 1.6490 1.6447 1.7419 1.6593
4.8024 4.8424 4.8026 4.8069 4.7783 4.7900
6.1793 6.2317 6.1804 6.1860 6.1492 6.1643
1,6300 1,6754 1,6301 1,6233 1o8636 1.6361
10 -4 0.01 1.00 0.01 0.05 0.i0 1.00
6.5209 6.5211 6.5227 6.5171 6.5164 6.5182 6.5374
3.0672 3.0584 3.1241 3.1338 3.0588 3.2313 4.2285
1.6490 1.6468 1.6783 2.3254 2.2505 2.3030 1.8368
4.8026 4.8052 4.7679 4.0503 4.1349 4.0967 5.0598
6.1805 6.1838 6.1358 5.2098 5.3197 5.2691 6.5336
1.6302 1,6372 1.6509 0.9726 1.0267 0,9943 1.5454
YI4*I02
YI5*I06
YI6*I0
YI7*I0
TOL
GEAR ii GEAR 21 EPISODE ii EPISODE 21 SIEU
34
Y13"I0-5
0.01 1.00 10 -4 0,01 1.00 1.00
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
5.0960 5.2375 5.0958 5.0748 5.8259 5.1146
8.7794 8.7990 8.7767 8.7815 8.6570 8.7923
1.3800 2.0011 1.3834 1.4799 1.9797 1.1932
4.4427 3.3321 4.439L 4.4043 4.0353 4.4675
10 -4 0.01 1.00 0.01 0.05 0.i0 1.00
1.9999 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
5.0964 5.1182 5.1610 3.0436 3.2123 3.1087 4.8445
8.7769 8.7755 8.8447 6.8001 6.9881 6.8661 8.1412
1.3840 1.4732 0.8152 1.7607 1.6760 1.8483 2.4357
4.4374 4.3936 4.5437 4.5987 4.3991 4.8281 7.7700
Environmental Software, 1988, Vol. 3, No. 1
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al COMPARISON OF QSSA AND GEARS METHOD FOR A REALISTIC AIR POLLUTION APPLICATION The analysis of air pollution chemistry problems, such as photochemical oxidant and acid rain formation, requires the solution of large sets of very stiff equations. For example, the analysis of 50-100 chemical species involved in one-to-several hundred reactions is not uncommon for these problems. Studies of chemical mechanisms typically make use of one of the Gear's methods(i.e., Eq. 2). However, when multidimensional analysis of the problems is performed, the computational burden of using Gear's method becomes too great, and one of the asymptotic methods is used. To better understand the trade-offs between the asymptotic methods and Gear's method a multiday box model calculation using the Lurman, Loyd, Atkinson [34] chemical mechanism (referred to hereafter in this paper as LLA) was performed. The LLA mechanism consists of 53 chemical species and 112 chemical reactions. For comparison purposes the Q6SA method as implemented in the STEM-II model (cf, Carmichael, et al., [35]) is compared with results obtained using the Gear's method. The QSSA method is as indicated by Eq. 18. However, the time steps are fixed at 30 seconds during the daytime , 60 seconds at night time, and 12 seconds at sunrise and sunset. The comparisons between Gear's method and the QSSA method are performed for the gas phase kinetics in a box model. In order to emphasize only the chemistry, the :~ physical processes such as emissions, depositions, dilutions, transports, etc. are not included. A test simulation adopted initial concentrations typical of polluted areas near sunrise: 40 ppb of NO, i0 ppb of NO 2, 20 ppb of S02, 1 ppb of 03, 30 ppb of alkanes, 8 ppb of ethene, 7 ppb of alkenes, I0 ppb of aromatics, 7 ppb of aldehydes and 1 ppb of ketones. Typical background concentrations of CH 4 and CO of 1650 and 125 ppb, respectively, were also included throughout this simulation. The ambient temperature and relative humidity were held to be fixed at 20 C and 50%, respectively. A multi-day simulation for the summer solstice at 40 N latitude was employed. The comparisons for key species are shown in Figure i. The percentage relative error based on results from the Gear's method is listed in Table 12. Except for relatively large discrepancies in species like PAN, the QSSA method predicts reliable results when compared with Gear's method. Figure 2 shows the percent relative error in the mass balance of sulfur and nitrogen. The use of the QSSA method often fails to maintain the mass balance. Moreover, the strong coupling between more reactive species result in numerical instability if a relatively long time step is used. Implementation of a linear transformation (or lumping) technique allows the circumvention of these two main difficulties. These procedures are explained in detail in Hesstvedt et ai.[33]. Results using the QSSA method with lumping are presented in Table 12 and Figures 1 ~ 2. Clearly, the lumping technique maintains the mass exactly for sulfur and within 0.5% for nitrogen throughout the simulation, except at the peak occuring around 18:00 of the first day. This seems to arise from numerical stiffness associated with a dramatic change at sunset. On the other hand, for the case without lumping, the strong discrepancies are seen in nitrogen compounds, which are..relatively more reactive than sulfur compounds. The QSSA method executes ~ l0 times faster than Gear's method for this application.
suited for specific types of applications. In order to assess and quantify the attributes of these techniques, five representative methods (i.e., GEAR, EPISODE, GERK, DE, and SIEU) were applied to test problems. These test problems were chosen to span the range of problem types, including: linear and nonlinear; small and large dimension; real and complex eignevalues; and slightly stiff and very stiff conditions. Results of the numerical test calculations indicate that: i.
If the accuracy of solutions is considered alone, GERK was the best candidate for nonstiff problems. Due to the small stability limit, this is not recommended to be used for stiff ordinary differential equations. DE was quite suitable for use in solving moderately stiff systems (SR < 102) with eigenvalues of Jacobian matrix close to the imaginary axis. For very stiff systems, SIEU was the best method at less stringent tolerance because of its short computation time. It worked better than GEAR and EPISODE if the dimension of problem is greater than ~ 15 - 20. However, it did not work as well as GEAR and EPISODE at more stringent tolerances. The two BDF methods, GEAR and EPISODE were roughly comparable in terms of accuracy and computing time, and each resulted in significant savings of time over the explicit methods DE and GERK. EPISODE resulted in the higher accuracy in each example and GEAR worked more efficiently, except in obtaining solutions under conditions when the eigenvalues were close to the imaginary axis. If it is necessary to choose one as an all-purpose stiff solver, EPISODE with routine 21 is recommended. For use in large air pollution models the asymptotic methods offer a good blend of accuracy and conservation properties of the solution.
2.
3.
4.
5. 6.
ACKNOWLEDGEMENT Partial support of this work was provided by California Air Resources Board through a subcontract with Sigma Research. REFERENCES [i]
Curtiss, C.F. and Hirschfelder, J.O. Integration of Stiff Equations, Proc. of National Adad. Sci., U.S.A., 1952, 38, 3, 235.
[2]
Dahlquist, G. G., "A Special Stability Problem for Linear Multistep Methods", BIT 3, pp. 27-43 (1963).
[3]
Eschenroeder, A. Q., Boyer, D. W., and Hall, J.G., "Nonequilibrium Expansions of Air with Coupled Chemical Reactions", Phys. Fluids, 5, NO. 5, pp. 615-624 (1962).
[4]
Amundson, N. R., "Some Further Observations on Tubular Reactor Stability", Can. J. Chem. Eng., 43, 2, 49-55 (1965).
[5]
[6]
Davison, E. J., "The Numerical Solution of Large Systems of Linear Differential Equations", AICh. E. J., 14, l, 46 (1968).
[7]
Mah, R.S.H., Michaelson, S. and Sargent, R.W.H., "Dynamic Behavior of Multi-component and Multistage Systems", Chem. ~ . Sci., I__~7,619-639 (1962).
[8]
Luss, D. and Amundson, N.R., "Stability of Batch Catalytic Fluidized Beds", A I C h . E . J . , 14, 211 (1968).
[9]
Slavickova, A., Hlavecek, V., and Kubicek, M., "Symposium on Computer in the Design and Erection of Chemical Plants", Karlory Vary (1975).
SUMMARY Many techniques are available for use in the numerical solution of stiff ordinary differential equations. However, each of the different techniques has its own relative strengths and weaknesses and is best
Amundson, N. R., Luss, D., "Qualitative and Quantitative Observations on the Tubular Reactor", Can. J. Chem. Eng., 46, 6, 425-433 (1968).
Environmental Software, 1988, Vol. 3, No. 1
35
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al [i0]
Willoughby, R.A., Ed., "Stiff Differential Systems", the IBM Research Symposia Series, Plenum Press, New York, NY (1974).
[ii]
[29]
Chan, Y. N. I., et al., "Solution of Stiff Differential Equations and the Use of Imbedding Techniques", Ind. En~. Chem. Fundam., 17, 3, 133 (1978).
Finlayson, B. A., "Packed Bed Reactor Analysis by Orthogonal Collocation", Chem. En~. Sci., 26, 7, 1081 (1971).
[30]
Finlayson, B.A., "Nonlinear Analysis in Chemical Engineering", McGraw-Hill, Inc. (1980).
(12]
Froment, G. F., "Fixed Bed Catalytic Reactors", Ind. En~. Chem., 59___, 2, 18 (1967).
[31]
[13]
Lee, J.C.M. and Luss, D., "The Effect of Lewis Numer on the Stability of a Catalytic Reaction", AICh. E. J., 16, 620 (1970).
Shieh, D. S-S, "The Evaluation of Numerical Techniques for Solution of Stiff Ordinary Differential Equations Arising from Chemical Kinetic Problems", M.S. Thesis, Chemical and Materials Engineering, University of Iowa (1983).
[32] [14]
Fleck, R.D., Kirwan, D.J. and Hall, K.R., "Mixed Resistance Diffusion Kinetics in Fixed-Bed Adsorption under Constant Pettern Conditions", Ind. En~. Chem. Fundam., 12, i, 95 (1973).
Young, T.R. and Boris, J.P., A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reaction flow problems, J. Phys. Chem., 1977, 81, 2424.
[33] [15]
Brayton, R.K., Gustavson, F.G. and Liniger, W., "A Numerical Analysis of the Transient Behavior of a Transistor Circuit", IBM J. Res. Develop., I0___, 4, 292 (1966).
[16]
Byrne, G.D. et al., "A Comparison of Two ODE codes: GEAR and EPISODE", Com~ and Chem. E., 1,133-147 (1977).
Hesstvedt, E., Hov, O., and Isaksen, I.S.A., Qusisteady-state approximations in air pollution modeling: comparison of two numerical schemes for oxidant prediction, Int. J. of Chem. Kinetics, 1978, x, 971. Lurmann, F.W., Lloyd, A.C. and Atkinson, R., A chemical mechanism for use in long-range transport/acid deposition computer modeling, JGR, 1986, 91, 10905.
[17]
Beinfeld, J.H., Lapidus, L. and Hwang, M., "Review of Numerical Integration Techniques for Stiff Ordinary Differential Equations", Ind. Eng. Chem, Fundam., 9, 2, 266 (1970).
[34]
[18]
Bjurel, G., Dahlquist, G.G., Lindberg, B., and Oden, S., "Survey of Stiff Ordinary Differential Equations", Report No. NA70-11, Department of Information Processing, The Royal Inst. of Tech., Stockholm, Sweden (1970).
[19]
Gelinas, R.J., "Stiff Systems of Kinetic Equations", J. of Comput. Physics, 9__, 222-236 (1972).
[20]
Enright, W. H., Hull, T. E. and Lindberg, B., "Comparing Numerical Methods for Stiff Systems of ODEs", Bit 15, pp. 10-48 (1975).
[21]
Lapidus, L. and Schiesser, W.E., Ed., "Numerical Methods for Differential Systems", Academic Press, Inc. New York (1976).
[22]
Warner, D. D., "The Numerical Solution of the Equations of Chemical Kinetics", J. Phys. Chem., 81, 25, 2329 (1977).
[23]
Fehlberg, E., "Low-order Classical Runge-Kutta Formulas with Stepsize Control and Their Application to Some Heat Transfer Problems" NASA TR R-315, July (1969).
[24]
Shampine, L. F. and Gordon, M. K., "Computer Solution of Ordinary Differential Equations", W. H. Freeman and Company (1975).
[25]
Preussner, P. R. and Brand, K. P., "Application of a Semi-implicit Euler Method to Mass Action Kinetics", Chem. En~. Sci., 36___,I0, 1633 (1981).
[26]
Gear, C. W., "Numerical Initial Value Problems in Ordinary Differential Equations", Prentice-Hall,F Englewood Cliffs, NJ, pp. 158-166 (1971).
[27]
Carmichael, G. R. and Peters, L. K., "A Regional Model for SO x Transport in the Troposphere", Atmospheric Environment, 18___,937 (1984).
[28]
Michelsen, M. L., "An Efficient General Purpose Method for the Integration of Stiff Ordinary Differential Equations", A I C h . E . J . , 22, 3, 594 {1976).
36
Environmental Software, 1988, Vol. 3, No. 1
[35]
Carmichael, G.R., Peters, L.K. and Kitada, T., A second generation model for regional-scale transport~chemistry~deposition, Atmos. Environ., 1986, 20, 173.
[36]
Gear, C.W., The automatic integration of ordinary differential equations, Cumm. Acm, 1971, 14, 185.
.~
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh eI al
=• .,.o
LEgENO
''iiiiiiiiiii
,,
L
°i
°0 _.e~
O * - MSIM 14/ LUMPING o - M S I M W/O LUMPINg ~.- gERR'S METHOD
./r
~f=
o 12
DRY
18
2t
6
12
18
24
ORT 2 TIME O" CLOCK
1
OR~( 3
i
12
IB
OR~
24
6
12
1B
24
ORY 2 I|ME O'CLOCK
1
ORY 3
y................ :I m @' SO= z
c~
c-I LEGEND + - MSIM W/ LUMPING o - MSIM W/O LUMPINg - gERR'S M~THOD
z
8
,
~
~ ,.=,@]
~ ~
LEGEND + - MSIM W/ LUMPINg o - MSIM W/O LUMPING
"
_,.4_=,~ ~,~,"; ;- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
¢0 ~_ .3 E. 6
,
,
,
,
,
,
,
12
18
24
6
12
18
24
DRY 1
OR'~ 2 T ] ME O' CLOCK
....
S
,2
DRY 3
,~
DRY
=,
1 TIME
Figure
I.
Comparison between Gear's method.
the Q S S A method
~
~
~
~'~
ORY 2 O' CLOCK
DRY
and
N ! TROGEN
z
el: _J ,0--
tr~
o o . H l TH LUHP ING
_
z
u~
w~
o
~
xx H I T H O U T
SULFUR
LUMPING
............. .,~
td >
......
J
.........
? e
,'2
,;
DRY l
Figure
2.
2'4
;
,'~
,'8
fi
DAY 2 T I'M[':O' CLOCK
Comparison of the mass balance with and without lumping.
s
DR1 3
-
Environmental Software, 1988, Vol. 3, No. 1
37
The Evaluation of Numerical Techniques: D.Shyan-Shu Shieh et al
Table I0 Numerical Results for Problem III w i t h u = 1.0 ....................................................................................................... TOL
GEAR ii GEAR 21 EPISODE 11 EPISODE GERK DE SIED
CPU TIME
LAST H
LAST ORDER
STEPS TAKEN
FCN EVAL
JACOB EVAL
1.0E-6
1.99
0.75
3
202
269
25
1.0E-6
1.57
0.83
5
159
191
24
1.0E-6
11.67
0.026
8
943
1510
144
1.0E-6 1.0E-6 1.0E-6 0.5 0.6
14.49 10.37 7.30 2.37 1.25
0.012 0.021 0.021 0.092 0.122
5 4(5) I0 1 1
1901 2787 927 907 289
2524 20798 1874 1814 578
227
Table II N u m e r i c a l S o l u t i o n Yi for Problem III at t=20 w i t h u = 1 ...................................................................................................... Y1
Y2
Y3
Y4*I09
Y5*I05
Y6*I0
...................................................................................................... GEAR Ii GEAR 21 EPISODE ii EPISODE 21 GERK DE SIEU (s=0.5) SIEU (s=0.6)
2.26E-8
1.46E-8
-2.311E-20
1.69
4.55
1.35
-1.13E-f1
-3.48E-II
2.16E-II
5.82
4.59
1.35
1.83E-87
-6.99E-88
1.805E-35
2.06
4.54
1.35
1.83E-87
-6.98E-88
1.805E-35
2.06
4.54
1.35
1.83E-87 1.83E-87 3.96E-66
-6.99E-88 -6.99E-88 -1.34E-66
1.805E-35 1.805E-35 2.683E-31
2.06 2.06 4.10
4.54 4.54 5.41
1.35 1.35 1.36
2.29E-60
6.82E-61
7.031E-30
5.40
5.83
1.37
1.26E-87
5.65E-88
1.850E-35
2.06
4.54
1.35
SOL'N
Table 12 C o m p a r i s o n b e t w e e n the QSSA Method and Gear's M e t h o d
Species
% Relative After 1 day MSIM(1) MSIM(2)
HNO 3 SO 2 Sulfate 03 ALKA ETHE AROM HCHO H202 PAN KET
0.24 -0.03 -0.62 -0.91 -0.24 0.49 -0.55 -8.22 -9.71 -11.5 -1.80
*
**
38
R e l a t i v e Error With Lumping
=
Gear
-
-5.96 0.ii -1.83 -2.50 -0.39 2.60 2.05 -9.45 -10.2 "~-15.6 -3.80
MSDIM)/Gear*I00
Environmental Software, 1988, Vol. 3, No. 1
Error* After 2 days MSIM(1)** (MSIM(2)*** 0.46 -0.16 0.45 -0.81 -0.59 0.39 -2.96 -6.80 -2.52 -13.3 -1.31
-5.70 0.i0 -0.75 2.57 0.53 4.79 1.73 -7.97 -3.67 -17.8 -3.30