Jerrold S.Rosenbaum
2382
Conservation Properties for Numerical Integration Methods for Systems of Differential Equations. 2 derrold S. Rosenbaum Department of Mathematical Sciences, Virginia Commonwealth University, Richmond, Virginia 23284 (Received May 5, 1977) Publication costs assisted by Virginia Commonwealth University
Consider an initial value problem for a system of s ordinary differential equations of the form y' = f ( t , y ) , y(0) = a where y = (yl,...,yJT, f ( t , y ) = (fi(t,y), ...,f& Y))~. If the physical system (and the corresponding modeling equations of the above form satisfies a linear constraint of the form &wiyi M where (zq) are constants (e.g., conservation of mass or number density) or a quadratic constraint of the form Cf=~C&,bjjyiyj= E where (bi,) are constants (e.g., conservation of energy), then it is desirable that the numerical method(s) used to solve the differential equation also conserve the same linear or quadratic constraint. Necessary and sufficient conditions for linear multistep methods and Runge-Kutta methods to conserve linear and quadratic invariants will be discussed. Also, if the integration method is implicit, the effect of the nonlinear equation solution method on the conservation of linear and quadratic invariants is discussed.
I. Preliminaries Consider a system of s differential equations
For (simple) mass action kinetics, the differential equations may take the form
dY l l d t = fl(t, Y 1, dY2W = f 2 ( t , Y 1,
d y j / d t = Pj(t, y ) - y j L j ( t , y )
Y 2 , * ' *, Y2, ' '
dYSldt = f,(t, Y 1 , Y 2 ,
*
'9
Y,) Y,)
.,Y,)
(1)
subject to the initial conditions Yl(t0)=
c1, Y 2 ( t O ) = c2, *
* ' 9 YS(t0)
= cs
(3)
y,L,(y), and S, are the production, loss, and other possible source terms, respectively. In one space dimension using 12 altitudes, 14 chemical species, and a cubic spline finite element discretization, one obtains a system of approximately 600 simultaneous ordinary differential equations. Ideally, if the system of differential equations, (3), satisfies one (or more) linear and/or quadratic conservation laws, (4)and (5), then the numerical solution which approximates the solution to eq 1 should satisfy the same linear and/or quadratic conservation laws. DEFINITION. Let (3) be a system of differential equations satisfying a linear conservation law (4) (respectively quadratic conservation law ( 5 ) ) , and let ((t",y")), n = 0, 1, 2, , .., be a discrete numerical solution to (3) produced by a given method. The numerical method is said to be linearly (respectively quadratically) conservative with respect to a given weight vector w (respectively with respect to a given symmetric positive definite matrix B ) if
...,Y , ) ~f, = (71, . . -,fJT,and co = (CI, .. .,
GIT.
A system of differential equations arising in chemical kinetics problems often satisfies one or more linear conservation laws (or linear constraints) and may satisfy one or more quadratic conservation laws (or quadratic constraints). The basic questions we shall address are as follows: (i) Suppose the system of differential equations, (3), satisfies one or more linear conservation laws of the form (4a)
or, by differentiating
0 (4b) for all initial conditions, where w = (wI, wz,.. ., w , ) is ~a vector of weights. Which numerical integration methods for eq 3 produce numerical solutions which satisfy the same linear conservation law&)? (ii) Suppose the system of differential equations, (31, satisfies one or more quadratic constraints of the form
WTf =
(5) for all initial conditions, where B = [b,] is a symmetric positive definite matrix. Which numerical integration methods for eq 3 produce numerical solutions which also satisfy the same quadratic constraints? The system of differential equations can arise directly as a system of ordinary differential equations or can be the result of applying a semidiscrete finite difference or finite element approximation to a system of partial differential equations (so called method of lines). yTBy = E
The Journal of Physical Chemistty, Vol. 81, No. 25, 1977
where P,(t,y ) is the production term and yiLi(t,y) the loss term for the ith speciesq4 An example of a partial differential equation is chemical kinetics transport in the atmosphere, where one has the equation1 a y i / a t = v * @ ( ~ Y~ i , i +) P ~ ( Y-) Y i L i ( Y ) + St (7) where vyi is the gradient of the species yi, and Pi(y),
where Y = 011, Y Z ,
wTy = M
(6)
(2)
In vector form, the above system may be written as dYldt = f ( t ,Y ) , Y(t0) = co
i = 1, . . .,s
...
wTyn = M n = 0,1,2, n = 0,1,2, (or (yn)TBy" = E
. . .)
(8)
(9) For linear conservation laws it can be shown that many numerical methods are linearly conservative (see section 11). This fact has several important ramifications. (i) The use of linear constraints as a check on the accuracy of the numerical integration method is often misleading. All you are generally checking for is round-off errors and not truncation errors (which are usually much larger). For example, if we let y" be the numerical approximation to y(t,) and f" = f(t,, y") then the "integration methods" yn+l
= Yn
(10)
Numerical Integration of Systems of Differential Equations
2383
and are both linearly conservative, though any relation between the answers they produce and the true solution is purely accidental. (ii) If a physical system satisfied a linear conservation law, then a numerical method which preserves the same linear conservation law should be used in order to eliminate one source of numerical error. The use of numerical methods which do not preserve linear conservation laws can lead to very poor and even wrong numerical solution^.^ (iii) One can choose to numerically integrate all s equations simultaneously or one can eliminate one (or more) differential equations by using the linear conservation laws. The choice depends upon, s, the number of equations in the system, and whether or not the elimination of one (or more) equations destroys the sparsity of the system (or the Jacobian). For small systems of differential equations (e.g., 10 or fewer equations) the elimination of one or more equations can generally be accomplished without encountering numerical difficulties. For medium and large systems of equations (e.g., 20 or more equations), it can cost more to eliminate one or more equations. The increase in work (or computer time) is proportional to the percent increase in the bandwidth of Jacobian and, in general, round-off error will also be increased. As a rule of thumb, if the Jacobian is dense, eliminating equations can be desirable, but if the Jacobian is sparse then eliminating equations will probably be undesirable. (iv) If the same numerical method is used to integrate all the differential equations then the numerical method will simultaneously conserve all the linear conservation laws. For quadratic conservation laws, it can be shown that only a few numerical methods are quadratically conservative (see section IV).
11. Numerical Integration Methods which Are Linearly Conservative Consider a general linear multistep m e t h ~ dof, ~degree
k
where yn-j = (yln-j,yzn-j, , , ., ~ p - j )f"-j~ , = f(t,-j, y"-j), and h is the time step. THEOREM l.ll If wTyL= M for 1 = h - I, n - 2, , ,, n k , and if (4)holds for all ( t ,y ) then wTyfl = M where y n has been computed by a consistent linear multistep method of degree k. Proof: See ref 11. A key item in the proof of theorem 1 was C$oaj = 0. This is equivalent to saying that the numerical integration method can integrate y' = 0 exactly. As such, we can extend theorem 1 to variable step linear multistep methods of the form8
.
where h, = t , - tn-l, the step size for the nth step, and (a'k-j] and { T k - j l are the coefficient for the nth step in the variable step method and are functions of {h,, hfl+ . .., hn-k+l)*
.
COROLLARY I. If wTyl = M for 1 = n - I, n - 2, . ., n k and if (4) holds for all (t,y ) then wTy" = M where y" has
been computed by a variable step linear multistep method
of degree K for which C$oaflj= 0 for all n. Thus, the variable step, linear multistep methods developed by Gear,5Hindmarsh and Byrne,' and Shampine and WattslGare all conservative. Hybrid methods8 of degree k of the form
can also be shown to be linearly conservative.ll Let us now consider a general rn-stage Runge-Kutta method5
yn = yn-1
+ mz
Ujhj
n
=
1,2, *
..
j= 1
(15)
where the kl satisfy the system of equations
k j = hf(t,-l
m
+ cjh, yn-' + 1=2
1
bjlh')j= 0. . ., m (16)
THEOREM 2. If wTyn-' = M and (4) holds for all ( x , y ) , then wTyn = M for n 1 1,where y" has been computed by an rn-stage Runge-Kutta method. Proof: If we multiply eq 15 by wT we obtain
~ T y =n W T y n - 1
+
m
j= 1
a.wTkj J
= ujryn-l
QED. Since Runge-Kutta methods are one step methods, all Runge-Kutta methods with variable step sizes are linearly conservative. Other methods which are linearly conservative include smoothing techniques (Lindbergg),averaging techniques (Liniger and Odeh"), and global and local extrapolation (Dahlquist,2 Gragg,6 Schryer13). The proof that these methods are conservative depends on two points: (i) all the numerical integration methods used to compute a solution point are conservative, and (ii) the linear combination of the underlying numerical integration methods used to produce the "final" numerical solution at each point is such that the sum of coefficients for the linear combination is 1.
111. Methods which Are Not Linearly Conservative Quasisteady state approximations4 are probably the largest class of numerical integration methods which are not linearly conservative. Depending upon which derivatives you set to zero, the remaining system of differential and algebraic equations may no longer possess the same (or any) linear conservation properties. A second method which is not conservative had been used by several authors (see ref 14 and 15). The numerical integration technique for eq 6 can be written as (yintl - y , " ) / h = P i n - yin+'lin
i = 1, . , ,,s (17)
where P," = P,(t,, y") and L? = L,(t,, y"). The problem with this method is that it evaluates the right-hand side of the differential equation a t two different time levels. For a more complete discussion of this method see ref 11. A third method which is not linearly conservative is the exponential Adams methods developed by the author of this paper.12 Consider the system of differential equations of the form
Y'=AY +&t,Y)
(18)
The Journal of Physical Chemlstty, Vol. 81, No. 25, 1977
Jerrold
2364
where A is an s X s matrix. One such exponential Adams method is yntl = e h A y n [hA]-l[ehA - Ilgn (19)
+
These methods were developed for differential equations describing a motion that is highly damped and for which g(t, y ) is discontinuous. The system of equatons does not have a linear conservation property.
IV. Methods Which Are Quadraticly Conservative Unlike the case for linearly conservative methods, there are relatively few methods which are quadraticly conservative. Suppose that for all initial conditions, the solution to eq 3 satisfies the quadratic conservation law V(y) = yTBy = E , a constant (20) where B is a positive definite symmetric matrix. An example of such a quadratic invariant can be energy in a closed nondissipative system. If we differentiate eq 20 we obtain
Since V ( y )is to be conserved, we require that for all y
V’(Y)f(Y) = 0 (22) For the special case V ( y ) = yTy = Cf,1yi2, then we obtain y T f = 0. This is equivalent to saying that the motion must be orthogonal (“perpendicular”) to the derivative of the change in the solution or that the solution must move on the surface of a sphere s dimensional space. THEOREM 3.3 Let V ( y ) be a real valued quadratic function. If every solution of the system of differential equations, (3), satisfies V(y(t))= E, a constant, then V W ) is constant also for every sequence produced by the implicit midpoint method (with either fixed or variable step). R o o t 3 The implicit midpoint method is defined by the equation y n t i = y n + hfntli2 (23a) where
(It should be noted that this method is A stable.) For a quadratic function, V ( y ) V(y””) - V(y”) = J;:+’V’(y) dy = V~(yntliZ)(yntlY”) (24) For the implicit midpoint method we, therefore, have V(yn+l)- ~ ( ~ =f lV’(yn+l/Z).h.fn+l/2 ) =0 (25) or V(ynt’) = V(y”) = constant QED. The above proof also yields the result that if V ( y )is a quadratic function and if V(y)is a nonincreasing function of t for every solution, then V’ ( y ) f ( y )I 0 for all y . The above proof relies on a special property of quadratic functions and on the midpoint rule and unfortunately does not lend itself to direct extensions. THEOREM 4. Let y ’ = Ay, where A is a skew symmetric matrix (AT = -A), and y T y = E , and let a numerical The Journal of Physical Chemistry, Vo/. 81, No. 25, 1977
S.Rosenbaum
method produce a solution ynfl = R(hA)y”,where R(hA) is a rational polynomial in hA. If [R(hA)]*R(hA) = I , the identity matrix, then the numerical method is quadratically conservative. Proof:
(yntl)Tyntl = ( Y ” ) [R(hA)ITR(hA)yn ~ = (y”)Tlyn n T n =(Y ) Y QED. Based on theorem 4,we now know that for linear systems, Oberechkoffs’ method and Butcher’s fully implicit Runge-Kutta methods that are m stage and order 2m are quadratically conservative. Unfortunately other types of Runge-Kutta rnethods or linear multistep method with degree 22 are not quadratically conservative. V. Solving Nonlinear Equations The equations that describe chemical reactions are generally stiff systems4s5of differential equations. That is, there is a wide variance in the time constants associated with the various chemical reactions. As such, an implicit numerical integration method designed for stiff systems should be used. Also it would be preferable if the iterative method, which is used to solve the nonlinear equations at each time step, would preserve the conservation property. DEFINITION. An iterative method for solving equations is said to preserve the linear conservation properly (respectively quadratic conservation property) if all iterates y b ) ,p = 0, 1, . . ., satisfy wTyb) = M where wT = (wl, . , ., LO,) (respectively ( ~ b ) ) ~ B = y bE)where B = [bill). It has been shown that for linear multistep methods and semiimplicit Runge-Kutta methods, the following methods for the solving nonlinear are linearly conservative: Newton’s methods, secant method, Broyden’s method, and sucessive substitution. For fully implicit Runge-Kutta methods, one must use a Gauss-Siedel-Newton (or Broyden or secant) or a Jacobi-Newton (or Broyden or secant) method in order that the nonlinear equation method preserve the linear conservation property.’l The corresponding results for the quadratic conservation property are not yet known to the author. VI. Concluding Remarks If the system of differential equations being solved satisfies a linear or a quadratic conservation property then a numerical integration technique which preserves the same quantity should be used. Lack of conservation can prove very embarassing (e.g., loss of mass). Also forcing conservation by “adjusting” the numerical solutions so as to satisfy the desired conservation property can lead to wrong answers. Finally, the fact that a particular numerical solution preserves the desired conservation properties is no indication whatsoever of the accuracy of the numerical solution (e.g., see the “methods” of eq 10 and 11). To calculate estimates of the errors in the solution other techniques should be ~ s e d . ~ p ~ , ~ J ~ References and Notes (1) J. S. Chang, A. C. Hindmarsh, and N. K. Madsen, “Simulation of Chemical Kinetics Transport in the Stratosphere”, in “Stiff Differentlal Systems”, R. A. Willoughby, Ed., Plenum Press, New York, N.Y., 1974. (2) G. Dahlquist, BIT, 3, 27-43 (1963). (3) G. Dahlquist and B. Llndberg, “On Some Implicit One-step Methods for Stiff Differential Equations”, Royal Instltute of Technology Report TRITA-NA-7302 (1973). (4) L. Farrow and D. Edelson, Int. J. Chem. Klnet., 6, 787-800 (1974). (5) C. W. Gear, “NumericalInitial Value Problems In Ordinary Differential Equations”, Prentice Hall, Englewood Cliffs, N.J., 1971. (6) W. B. Gragg, S I A M J Numer. Anal., 2 , 384-403 (1965).
Analysis of Multiparameter Model Systems (7) A. C. Hlndmarsh and G. D. Byrne, "A Polyalgorithm for the Numerical Solution of Ordinary Differential Equations", Lawrence Livermore Laboratory Report UCRL-75652, 1974. (8) J. Lambert, "ComputationalMethods in Ordinary Differential Equations", Wlley, New York, N.Y., 1973. (9) B. Lindberg, BIT, 11 29-52 (1971). (10) W. Liniger and F. Odeh, IBMJ. Res. Dev., 16, 335-348 (1972). (11) J. S. Rosenbaum, J . Comput. Phys., 20 259-267 (1900). (12) J. S. Rosenbaum, "Implicit and Explicit Exponential Adams Methods for Real Time Simulation", IMACS Symposium on Simulation Software arid Numerical Solution of Differential Equations, 1977, to be published.
2365 (13) N. L. Schryer, "A Users Guide to WDES, a Double Precision Ordinary Differential Equatian Solver", Bell Laboratories, Computer Science Technical Report No. 33, 1975. (14) T. Shimazakl and A. R. Laird, J . Geophys. Res., 75, 3221 (1970). (15) W. Stewart and M. J. Hoffert, "Stratospheric Contamination Experiments with a One-Dimensional Atmospheric Model", International Conference on Environmental Impact of Aerospace Operations In the High Atmosphere, Denver, Colo., June, 1973. (16) L. Shampine and M. Gordon, "Computer Solution of Ordinary Differential Equations-The Inkial Value Problem", W. H. Freeman, San Francisco, Calif., 1975.
Nonlinear Sensitivity Analysis of Multiparameter Model Systems R. I. Cukler,+ Department of Chemistry, Michigan State University, East Lansing, Michigan 48824
H.
B. Levlne,
Jaycor, 140 1 Camino Dei Mar, Del Mar, California 920 14
and K. E. Shuler" Department of Chemistry, University of California- San Diego, La Jolla, California 92093 (Received April 26, 1977) Publication costs assisted by the National Science Foundation
Large sets of coupled, nonlinear equations arise in a number of disciplines in connection with computer based models of physical, social, and economic processes. Solutions for such large systems of equations must be effected by means of digital computers using appropriately designed codes. This paper addresses itself to the critically important problem of how sensitive the solutions are to variations of, or inherent uncertainties in the parameters of the equation set. The sensitivity analysis presented here is nonlinear and thus permits one to study the effects of large deviations from the nominal parameter values. In addition, since all parameters are varied simultaneously, one can explore regions of parameter space where several parameters deviate simultaneously from their nominal values.
Synopsis Sets of coupled, nonlinear equations arise in a number of disciplines in connection with computer based models of physical, social, and economic processes. These sets of equations may be differential, integral, or algebraic. They arise in such widely different fields as reaction kinetics, combustion, air pollution, weather forecasting, upper atmosphere phenomena, seismic analysis, operations research, systems analysis, economics, etc. These model systems may contain as many as 100 equations, a very large number of parameters (in the form of coupling coefficients such as rate coefficients, transport coefficients, economic coefficients, etc.) and a very large number of output variables. Solutions for such large systems of equations must be effected by means of digital computers using appropriately designed codes. As is well known, the computer solution of such large sets of equations can be quite expensive. Even after such solutions have been achieved one is still faced with a critically important problem: How sensitive is the solution t o variations of, or inherent uncertainties in, the p a rameters of the equation set? This problem of sensitivity 'Alfred P. Sloan Fellow. Supported in part by NSF Grant No. CHE 74 09442. *Supported in part by the Advanced Projects Agency of the Department of Defense, monitored by the U. S. Office of Naval Research under Contract No. NOOO14-69A-0200-6018during the period Aug 1,1972-July 31, 1974.
is central to the understanding of the behavior of systems, and of the models representing such systems, which contain a large set of coupled equations. It is clearly important to know how sensitive the output variables are to changes of, or uncertainties in, the parameters and which of the variables are sensitive (or not sensitive) to which of the parameters. Until this information is available, any proposed model must be suspect as a valid representation of the real system. Furthermore, the accuracy to which the model parameters, Le., coupling coefficients, need to be determined via calculation or experiment depends upon the sensitivity of the output variables to the value of the parameters. In addition, any desired optimization of various output variables with respect to the coupling coefficients requires a knowledge of the sensitivity. Any attempt to determine the sensitivity by solving the set of equations over and over again, varying one parameter at a time over a series of values while holding all the other parameters fixed a t some specific values, becomes prohibitive in time and expense for the large systems discussed here. This is readily demonstrated by a simple calculation. For a model system of many coupled differential equations with n parameters and m output variables, the above procedure for z different values for each of the parameters would require zn integrations for each of the m variables, i.e., a total of m(z)" integrations. For m, z , and n large, not only will the computations be prohibitively expensive, but the printouts would be so numerous that the analysis The Journal of Physical Chemistry, Vol. 81, No. 26, 1977
