Application of Backward Differentiation to Chemical Kinetics
2409
Extrapolant Formulation of the Backward Differentiation Method with Application to Chemical Kinetics Equationst A. M. Wlnslow Lawrence Livermore Laboratory, University of California, Livermore, California 94550 (Received May 27, 1977) Publication costs assisted by Lawrence Livermore Laboratory, University of California
The implicit backward differentiation method of Curtiss and Hirschfelder as generalized by Gear to systems of ordinary differential equations takes its simplest form when extrapolants to the solution at each step are used for the history array. The extrapolants, which form successive approximations (of order 0 to k 5 6) to the solution at each step, are easily obtained recursively and are used for prediction, correction, interpolation, and step and order control. Change of the time step at each step is assumed. Examples are given of applications of this method to systems of chemical kinetics equations.
I. Introduction
which permit use of a time step At
The backward differentiation method was invented in 1952 by Curtiss and Hirschfelder' for the numerical integration of "stiff" ordinary differential equations, e.g., chemical kinetics equations, for which the older numerical methods required impractically small time steps. The method was extended to systems of equations by Gear' in 1968, who developed a practical scheme with automatically varying order and time step. Since that time other formulations of the method have been developed by Brayton Byrne and Hindet al.,3 wins lo^,^ Rubner-Peter~en,~ marsh,e Van Bokhoven? and Krogh.8 Although essentially equivalent to each other, these formulations differ in the manner in which they store the "history array", Le., the previous solution values, which leads to quite different algorithms. The extrapolant formulation presented here4 is simpler than any other formulation and an interested scientist can easily write his own computer program based on it. This formulation was apparently rediscovered by Van Bokh ~ v e nand , ~ a somewhat similar method has been developed by Rubner-Peter~en.~ In this paper a summary of the method is given together with some examples of applications to chemical kinetics. Derivations and further details may be found in ref 4 and in a forthcoming paper.g
111. Difference Equations
11. Differential Equations We consider a first-order ordinary differential equation
(1)
dY / d t = f(Y
or a system of such equations dyi/dt = fi(y') i = 1,2, . . . , I (2) where y' represents the set of variables yl, y2, ..., YI. In chemical kinetics system (2) usually can be put in the form dyi/dt =-ai(y')yi
+ bi@)
i = I, 2, . . .,I
Work performed under the auspices of the U.S. Energy Research and Development Administration under Contract No. W-7405-Eng-48.
or aiAt
>> 1.
Let us temporarily restrict our discussion to eq 1,and assume that we have calculated the solution yo, y', ...,y" at times (not necessarily equally spaced) to, t', ..., t", and we wish to calculate yn+' by the backward differentiation method. We begin by constructing an interpolating polynomial y p ( t )of order h which passes through the unknown value ynt' together with the k past values y", y"-l, ...,ynMk+' The backward differentiation method consists simply of using the derivative of y,(t) at tn+'for the left-hand side of eq 1:
.
dy,/dtl,=,n+l
= f(y"+'
1
(4)
This usually leads to an implicit equation for yntl which is solved by an iteration method. The most convenient form for y,(t) is the Newton interpolating polynomial which contains the past values (y", ) the form of divided differences. As shown y=-1,...,Y ~ - ~ + 'in below, there are several advantages to rewriting these divided differences in terms of the extrapolants (yont1, ylnt', ..., yk-ln+l)which contain exactly the same information in a different form. Each extrapolant yj"+' is defined to be the value of yntl obtained by extrapolating a polynomial of degree j from past values. Together they form a sequence of successive approximations to the unknown yntl. The lowest order extrapolant is yo"+l = Yn (5) and the others can be easily obtained from the recursion relation4 Yk+1
n+l =
ykn+' + Tk(yn - ykn)
(6)
where y n is the calculated value of y at t", ykn are the previous extrapolants a t tn,and
(3)
where bi and -aiyi represents production and loss terms, respectively, for chemical species y'. The quantity acl is the characteristic loss time for yi. Stiffness exists when, for at least one yi,the characteristic time over which y Lis changing yi/yi >> a;'. In this case we use stiff methods
>> u;',
(7) When the Newton interpolating polynomial is written in terms of the extrapolants and then differentiated, eq 4 takes on the simple form4 kgl(:+l j=o
n+l - n-j
= f(yn+i)
(8)
The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
A. M. Winslow
2410
This is the basic equation of our method. For k = 1, eq 8 becomes (yn+' - y n ) / A t = f ( y n + ' )
(9)
where At = tn+l- tn, which is known as the backwards (or implicit) Euler method. For k > 1,the succeeding terms on the left-hand side of eq 8 represent successively higher order corrections which increase the accuracy. Applying eq 8 to system (3) and solving for yinfl we obtain
3
L
4
u3
X
I
h- 1 L
rjYi,jn+l + bin+lAt
I.
yin+i = j = O
12-1
C rj + ain+'At
(10) llt
j =1
i = 1, 2,. . .,I
c
where yiJn+' is the j t h order extrapolant to ylntl and tnt1 - n
rj = t n + l
t
-10
0
Log t ( s e c )
- tn-j
Figure 1. Methane oxidation: [CH,]
We have written a;+l and bln+lto indicate that they are functions of the yln+l. The stability of the backward differentiation method for chemical kinetics equations is suggested (although not proved) by eq 10. In the limit of large time step y?+lbi/a,, the quasisteady solution which thus serves as a bound. Furthermore if a, > 0 and b, > 0, as they usually are, we see that ylntl > 0 provided that ylJn+l> 0. Because of the separation of loss and production terms, eq 10 are in a form suitable for solution by functional iteration. From initial guesses for the yln+lwe calculate a, and b,, then recalculate yln+land continue until convergence is reached. This method has the advantage over iteration by Newton's method that no derivatives or Jacobians are needed, and no linear systems need be solved, and it is somewhat faster than Newton's method. It has been used successfully for large (30 species) stiff systems. However its convergence is problem dependent, although various devices have been found which extend the domain of convergence. For very large systems or for reactive flow calculations, computer storage and speed limitations make it advantageous to avoid Jacobian calculations, and more work needs to be done to improve the convergence of this functional iteration method. Equations 10 are also suitable for iterative solution by Newton's method. Since they are in "row-scaled'' form, the main diagonal of the Jacobian matrix is all ones, which helps to improve the conditioning of the matrix. The summations are constants in the iteration, and the Jacobian makes use of the derivatives aa,/ay, and ab,/ay,.
IV. E r r o r Control This formulation of the backward differentiation method uses a continuously variable time step and order, based on control of the local truncation error. The leading term in the local truncation error in dy/dt is just the first neglected term on the left-hand side of eq 8, namely, ( y n + l - ykn+l)/(tn+'- Pk), and the local truncation error in yn+' can be shown to be essentially the product of this with At. In practice we use simply yn+l- ykn+', so that the relative error at each step is taken to be (in Euclidean norm)
It can be shown that
0 ' i - L L - L -30 -2 0
6k
= (constant) X
The Journal of Physical Chemistty, Vol. 8 1, No. 25, 1977
Thus at
vs. log time.
each step the time step Atnewfor the next step is calculated from
Atnew = (60/6h)1i(ht1),it (12) where 60 is the desired relative local truncation error, which I6o 5 is usually in the range By computing Atnewfor k - 1and k 1,and comparing these values of Atnew,we have a criterion for changing the order h by fl whenever it will lead to a larger step. Changing the order merely means changing the number of terms in the sum on the left-hand side of eq 8 or in the numerator and denominator of eq 10. As with other implementations of the backward differentiation method, various controls on the change of k and At, somewhat problem dependent, are included in a practical code. For example, if bk exceeds bo by more than about 3070,the entire step is repeated with a smaller value of At.
+
V. Advantages of t h e Method The main advantage of this formulation is the economy obtained by using the extrapolant history array. The extrapolants are used as predictors (eq 6) and they also carry all the information needed for correction (eq 81, control (eq ll),and interpolation (for output purposes). If this formulation is compared with others, the gain in simplicity is striking.
VI. Numerical Results A computer code (named KIN) has been written for the solution of systems of chemical kinetics equations of the form (3) by the method described above. It has been tested on several chemical systems and the results compared with calculations on the same systems with the programs GEAR'' and EPISODE.^ Some results are shown in Figures 1-11 which were produced by the computer, each dot representing one time step. The problems are as follows: (1) Reaction of a stoichiometric mixture of methane and air at 1 atm and 2000 K, using the reaction mechanism of Bowman, Pratt, and Crowell which has 16 species and 24 reactions. (2) A one-equation simulation of diurnal photochemistry.12 (3) The Belousov reaction mechanism of Field and Noyes,13 which has three species.
Application of Backward Differentiation to Chemical Kinetics
81
2411
"I
I
~
2
1
1
-3 0
(sec)
Log t
Flgure 2. Methane oxidation: [CHO] vs. log time.
,
, I , ,
, ,
,
, , I ,
-2 0
-1 0 Log t ( s e c )
,
I~
0
Figure 5. Methane oxldation: order vs. log tlme.
1
I I
1.0
t
I v
I
n
I
UY N 4 0
N
0
z
:-l6I L
0
,
0.5
c
-18
~
-2 0
L l
-2 0
-30
-1 0
Log t
0.0
0
2.0
1.0
3.0
t ( l o 5 sec)
(sec)
Flgure 6. Diurnal photochemistry simulation: concentration vs. time.
Figure 3. Methane oxidation: log [NO2] vs. log time.
'11
U
v1
v
c,
c,
a
-I
-2
-30
t
1 LI
~
-30
u
i
,
1
.
,
1
~ 1
~1
1~
, , ,
~ ~
,
,,,[I
-1 0 Log t (sec)
-20
Figure 4. Methane oxidation: log time step vs. log time.
0
1
0.0
-
1.0
t (lo5
L
2.0
3.0
sec)
Figure 7. Diurnal photochemistry simulation: log time step vs. time. The Journal of Physlcal Chemistry, Vol. 81, No. 25, 1977
A. M. Winslow
2412
t
5;
d
L
E -4k
”
t
1r
0 -I0
I
I
1
2
,
_
.
L
I
A
3
5
0
10
t (sec)
t ( l o 5 sec) Figure 8. Diurnal photochemistry simulation: order vs. time.
Figure 10. Belousov reaction: log time step vs. time.
I
5-
t
4L L a
-0
+ ,
,
& 3-
-
1
I _ _ _ L _ I L _ _ _ L l - . 2
0
5 t (sec)
10
C
5
10
t (sec)
Figure 9. Belousov reaction: log [HBr-] -t 6.52 vs. time.
Figure 11. Belousov reaction: order vs. time.
The three codes gave results for these particular calculations that agreed to a t least four figures, with computing speeds that differed by less than 20%. These problems were chosen to illustrate the stability and efficiency of the method, which are reflected in the size of At, the amount of change of At and k , and the number of steps that needed to be repeated. In these respects each problem presented a different challenge: Problem 1 has several species concentrations of greatly different magnitudes and rates of change. The initial time s s) was reduced by the program to At = step for accuracy. The time step then increased steadily through the induction period, by a factor of 10 at each step while k = 1,and by smaller factors for k > 1. At t = s ignition took place, with rapid changes of most of the species (see Figures 1-3). The timestep was reduced somewhat (Figure 4) but then resumed its increase up to t = 1s. The order increased in a step-like fashion (Figure 5 ) until t = 0.1 s when roundoff due to the attainment of equilibrium (right-hand sides of eq 3 equal to zero to eight figures) introduced noise which limited the further increase of the timestep and caused the order to be reduced to unity.
Problem 2 has a periodic solution that is nearly a square wave with a 24-h period (Figure 6) which is difficult to integrate accurately because of the sudden large time step changes at the “corners” of the square wave. The time step increased to the maximum allowed (At,, = 2 h) during the flat portions of the cycle (see Figure 7). In order to reproduce the rapid changes at sunrise and sunset, which took place over a few milliseconds, the timestep was automatically reduced by about five orders of magnitude, requiring five or six repetitions of the step since the extrapolated values were no longer accurate. The order also exhibited periodic changes (Figure 8) but no instabilities occurred. The GEAR program, however, when run with certain parameter combinations, was unable to complete this pr0b1em.l~ Problem 3 also exhibited sudden time step changes combined with nonlinear interactions of the reacting species (Figures 9-11). No instabilities occurred and no repetition of steps was necessary.
The Journal of Physical Chemistry, Vol, 8 1, No. 25, 1977
VII. Conclusion A new formulation of the backward differentiation method has been described which is equivalent to, but
Chemical Kinetics Modeling Based on the Taylor Series Expansion
simpler than, other formulations. Because of its simplicity it should be useful to those who wish to write their own chemical kinetics programs or to incorporate chemical kinetics into a reactive flow calculation. The basic equations are eq 6, 10, 11, and 12, together with a subroutine to calculate the aiand bi and their derivatives, and a linear equation solver to carry out the iterative solution of eq 10. For computers with small word length, roundoff error can be reduced by replacing yn+land ykn+' by the backward differences Aynfl = yn+' - yn and Aykn+' = ykn+' - yn. In this case eq 6 becomes
so that the extrapolation is done with backward differences. Then the extrapolants are recovered from the relations ykn+' = yn + Aykn+'and the rest of the calculation proceeds unchanged. References and Notes (1) C. F. Curtiss and J. 0. Hirschfelder, Proc. Natl. Acad. Sci., 36,735 (1952). (2) C. W. Gear, Inf. Process., 68 (1969).
2413 (3) R. K. Brayton, F. G. Gustavson, and G. D. Hachtel, Proc. I€€€, 60 (l), 98 (1972). (4) A. M. Winslow, "Numerical Integration of Chemical Rate Equations", Lawrence Livermore Laboratory Report UCRL-73911 (Jan 1972) and "A New Formubtionof Backward Differenced Linear Mutisteo Methods for Ordinary Differential Equations", Lawrence Livermore Laboratory Report UCIR-681 (Jan 1973). (5) T. Rubner-Petersen. "An Efficient Alaorithm Usina Backward Time-Scaled Differences for Solving stiff Differencal-Algebraic Systems", presented at the Seminar in Numerical Analysis, Royal Institute of Technology, Stockholm (Oct 10-12, 1973). (6) G. D. Byrne and A. C. Hindmarsh, ACM Trans. Math. Software, 1, 71 (1975). (7) W. M. G. Van Bokhoven, I€€€ Trans. Circuits Syst., 22, No. 2, 109 (1975). (8) F. T. Krogh, "Changing Stepsize in the Integration of Differential Equations using Modified Divided Differences", in Proceedings of the Conference on the Numerical Solution of Ordinary Differential Equations, University of Texas at Austin (Oct 1972). (9) To be submitted to the J. Comput. Phys. (IO) A. C. Hindmarsh, "GEAR: Ordinary Differential Equation System Solver", Lawrence Livermore Laboratory Report UCID-30001, Rev. 3 (Dec 1974). (1 1) B. R. Bowman, D. T. Pratt, and C. T. Crowe, Symp. ( I t . ) Combust., [Proc.], 14th, 824 (1973). (12) G. D. Byrne, A. C. Hindmarsh, K. R. Jackson, and H. G. Brown, "A Comparison of Two ODE Codes: GEAR and EPISODE", University of California, Lawrence Livermore Laboratory Report Preprint UCRL-79141, January 1977; to appear in Comput. Chem. Eng. (13) R. J. Field and R. M. Noyes, J . Chem. Phys., 60, 1877 (1974).
A Numerical Method for Chemical Kinetics Modeling Based on the Taylor Series Expansion John
P. Kennealy and William M.
Moore*
Department of Chemistry and Eiochemistry, Utah State University, Logan, Utah 84322 and the Air Force Geophysics Laboratories (OPR), Bedford, Massachusetts 0 1730 (Received May 27, 1977) Publication costs assisted by the Air Force Geophysics Laboratory
The sets of non-linear, first-order differential equations which describe the time behavior of chemical systems are ordinarily integrated with numerical techniques that utilize the analytic form of only the first derivative of species concentrations. However, these differential equations are somewhat unique in that the analytic form of all higher order derivatives can also be generated quite readily. A numerical integration algorithm based on the Taylor series expansion has been developed which exploits the specific form and symmetries of the differential equations governing chemical kinetics to provide an ease and efficiency of operation not found with other techniques. The algorithm is self-starting and a local convergence test is applied to vary the step size through the integration. In addition, a stiff equation technique is folded into the solution in the form of perturbed derivatives for any component which, by transient pseudo-steady-state behavior, forces the integration to small step sizes inconsistent with the time constants of the reaction set.
Introduction Over the last decade or two the increased availability of good kinetic data and the rapid improvements in digital computer technology have made the numerical modeling of large scale.chemica1 systems a most fertile area of activity. For example, numerous models have been developed and are being applied in the study of atmospheric chemistry problems, ranging from local air pollution effects through stratospheric ozone depletion to ionospheric perturbations produced by aurora. During the same time period, many other disciplines have developed their own requirements for sophisticated numerical modeling of a broad spectrum of phenomenology and descriptive systems. Much of the spectrum of modeling requirements shares
* Address correspondence t o this author at U t a h State University.
the same general concern, i.e., the numerical integration of sets of coupled differential equations. Consequently, these apparently common requirements have been a major incentive for a large body of work directed toward the development of better general numerical integration techniques. Today one sees a profusion of elegant techniques being applied to the large variety of modeling problems. Whereas not long ago the chemist concerned himself with a relatively few extrapolative and predictor-corrector methods (e.g., Eluer,' Runge-Kutta, Kutta-Merson,2 etc.), today he sees numerical modeling being increasingly dominated by rather complex, general purpose methods such as the Gear a l g ~ r i t h m . ~ , ~ These multistep implicit numerical integration methods, which usually achieve significant reductions in execution time for typical "benchmark" problems, require the construction, handling, and inversion of matrices; consequently, storage requirements for the matrix methods The Journal of Physical Chemistry, Voi. 81, No. 25, 1977
