255
Ind, Eng. Chem. Fundam. 1981, 20, 255-266
Adaptive Semi-Implicit Runge-Kutta Method for Solution of Stiff Ordinary Differential Equations George J. Prokopakls and Warren D. Selder” Department of Chemical and Biochemical Engineering, Universb of Pennsylvania, philade@hla, Pennsylvania 19 104
Parameters of a second-order, noniterative, A-stable, semi-implicit Runge-Kutta method are adjusted to improve 0s accuracy. The step-size, h, is projected and a pseudo-eigenvalue is estimated for the linear ODE, dy,ldt = &ye,using values of the stiff variable, ya over two previous time-steps. Here, y8 is the variable having the largest fractional rate of change. The product M, is used to adjust the parameters of a semi-implicit Runge-Kutta method to ghre improved accuracy in tracking the stiff variable. The adaptive a l g o ” has been demonstrated to Integrate accurately linear and nonlinear stiff systems with and wbout oscillations. After limlted testlng, the algorkhm appears to be competitive and in some cases more efflclent than higherorder, seml-impllcit Runge-Kutta methcds and the GEAR program. Three of the test systems of ODES represent chemical reaction systems including the Belousov limit cycle.
Introduction The initial value problem
and the solution Y1
= YlwX’t
y2 = CleXlt+ CzeXlt arises often in analysis of chemical reaction systems and unsteady distillation towers, among numerous engineering applications. Usually, the m ordinary differential equations (ODES) are nonlinear and eigenvalues of the linear approximation dY dt z JY
(2)
range over several orders of magnitude. J is evaluated at times ti and the approximate solution about ti is y1 z dlleXlt+ ... + dlmeXm*
ym z dmleXlt+... + dmmeAmt
t 2 ti
(3)
...
where X1, ..., A, are eigenvalues and the 4 = [ d l j dmjIT, j = 1, ...,m, are eigenvectors. Most physical processes dampen disturbances and their models (ODES) respond by decaying to a steady-state solution. The ODEs are commonly referred to as “inherently” stable. For these systems, with complex eigenvalues Xi = XRi f iXIi
(4)
A, = RefA,:]< 0, j = 1, ..., m, and AIj is the frequency of oscillation in the transcendental terms of the solution; sin AIjt and cos AIjt. For XIj = 0, no oscillations occur as the solution decays to a steady state. Stiff ODEs. Aiken and Lapidus (1974)present two linear ODES with real eigenvalues
0198-4313/81/1020-0255$01.25/0
(64 (6b)
They plot the solution for yz with IXzl 2.732. Therefore, a Rosenbrock formula, with different coefficients, having third-order accuracy, cannot be strongly A-stable, but strong A-stability is desirable when hlRe(X)J, becomes large in stiff systems. To achieve an algorithm that approaches strong A-stability as h(Re(XJI,, increases, we begin by setting a = al = a2 and c1 = 0 to obtain a second-order approximation with one Jacobian evaluation and one LU factorization per time-step. The parameters are related by w1 + w2 = 1 first order (41) a + w2b = 1 / 2 second order (42) 1
Y=
+ (1- 2a)hX + (a2- u + ~ 2 b ) ( h X ) ~
(43)
(1- ahX)2
+~
a2 - u Ym
=
2 b
a2
(44)
We have three equations (41,42, and 44) and five variables (w,, w2, a, b, y-); b is related to w2 using eq 19 w2b2 = 1/3
(45)
which matches third-order terms, in the hope of improving the second-order approximation. One degree of freedom remains and we select ym,a parameter we adjust to approach zero as hlRe(X}lbecomes large. We relate ymto a by substituting for w2 in eq 42 to give b=
1/3 -a
and substituting in eq 44 a2
[ = r - 1, r, r
ym =
+1
where the parameters ml,m2,and m3 are a function of E. In addition, [ at the minimum is taken as the order of the multistep formula for the next step. Normally, only one iteration of a Newton-Raphson method is needed to solve for K + in ~ eq 34 and J is reevaluated once every few time-steps. Adaptive Algorithm Like the multistep methods, we have developed a strategy for adapting an integration formula to give improved accuracy, but we use a low-order, noniterative, A-stable, semi-implicit Runge-Kutta method. First, we introduce the method before presenting the adaptive strategy. Extension of Rosenbrock Algorithm. Returning to the Rosenbrock algorithm, eq 13-15, the limit of characteristic root, eq 20, as hlRe(X)I m is ala2- wla2 - waul + w2bl ym = (36)
-
ala2
Combining eq 16, 17, and 36 a, + a2 - a,a,(l
- Ym)
= 1/2
and combining eq 17 and 18 a, + a2 - 2ala2 = 1/3 Eliminating a, in eq 37 and 38 6(1 + Y ~ ) -u2(2 ~ +~ ym)a2+ 1 = 0 Hence, we see that for real a2
(37) (38) (39)
- 2a
+ 1/2
a2
(47)
In summary, parameters for the second-order approximation are expressed in terms of ym 1 a=
+ [l- (1- ym)/21”2 1- y m b=
~2
1/3 -a
= 3((1/2) -
(48)
w1 = 1 - wp
For example, when ym= -0.25, a = 1.29, b = -0.422, and x+l= J’i - O.872k1 + 1.872k2 (49) and when ym= -0.20, a = 1.36, b = -0.388, and x+l= E - l.216k1 + 2.216k2 (50) Local Truncation Error. To estimate local truncation error, we introduce an imbedded partner y’i+l = yi + ~ ’ 1 k + 1 ~’2k2 (51) with coefficients matched for first-order accuracy W‘l+ w’2 = 1 (52) a and b are unchanged so that kl and k2need not be recomputed and a2 - u + whb 7” = (53) a2
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981
250
I.
-Y a ym=-O.l Y,
s
- 0.25
Figure 3. Characteristic root for second-order formulas as a function of -hX, and parameter y-.
/'
,
-I.
Solution:
I
Figure 2. Linear approximation for the stiff variable, ye, at ti.
To give good estimates for large hlRe(X)I,ym'is equated to zero, and a - a2 WlZ = (54) b That is, the imbedded partner is strongly A-stable. Subtracting eq 51 from eq 13 we obtain an estimate of local truncation error 4 = Yi+i - &+I = (Wi- W'i)ki + (Wz- W'2)kz (55) = (wz- wh)(kz - 4) w(k, - k1) where
w = 3((1/2) - a ) ( d - 2a + (1/2)) (56) Tracking the Stiff Variable. Before completing a time-step, we select the stiff variable having the largest I(.~j,i+~ - yjj)lyj,il, j = 1, ..., m, and compute the pseudoeigenvalue, &, for the linear approximation dY, -= (57) dt the solution of which is illustrated in Figure 2. It remains to show how to compute but first we justify use of this approximation. Consider integration of eq 57 using the Rosenbrock formulas with parameters in eq 48. For each choice of ym, we obtain a poor approximation to the solution of eq 57 except over a small range of -h&. This is shown in Figure 3 where the characteristic root for each formula is plotted as a function of -h&. Since the exact solution of eq 57 is ys(t.r + l 1 = y8(t.)eL(t,+1-ti) I = ys(ti)ehL (58) and the solution of each difference equation is Ysj+l = 7(hUYSj (59) when the difference equation has y(h&) = ehL (60)
a,
Dl 0
2
6
L
0
- ht,
10
12
11
18
Figure 4. Locus of points where eq 13-15 with pcamebre in eq 48 give exact solution for the linear ODE, dy, f d t = 0,.
it gives the exact solution. Hence, given h& (it remGns to estimate As), ymcan be selected to give the difference equation that intersects the curve traced by eq 60. This locus of intersections is plotted in Figure 4 along with the empirical equation A +E ym = (61) B (hQ2 D + h&
+
+-
where A = -4.9221, B = 21.1642, C = 0.5287, D = -0.6889, E = 0. Note that these parameters give an excellent approximation for -hX > 3, the region of importance in integration of stiff systems. For integration of eq 57, accuracy of the numerical solution is limited by eq 61; the exact solution corresponds to a perfect fit. Adapting the Integration Formula. With an estimate of X,, ymcan be computed for accurate integration of eq . approach is to use this X, and an 57 to give ~ , , i + ~ Our estimate of stepsize, h, to compute ymand the parameters (eq 48) for integration of the 5ystem of ODESover the next step. We obtain improved accuracy in tracking the stiff variable, as well as the remaining variables. Estimating the Pseudo-Eigenvalue, A,. Next, we show how to estimate the pseudo-eigenvalue for the stiff variable. The ODE for the stiff variable is
and its linear approximation is
?!dt !
- f ( ")yj ayj N
j=1
m
= j=lJajyj
(63)
where JSj is an element of the row s and column j in the
260
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981
Jacobian matrix. Combining the definition for the pseudo-eigenvalue dY, P xy, (57) dt with eq 63 and differentiating with respect to ya
When the derivative is estimated using first-order approximations
Substituting in eq 64 and multiplying both sides by h t,.2t,.l
t,
17.1
t
The vector kl is given by [I - haJ]kl = hf(yi]
Figure 5. Extrapolation for step-size with linear and quadratic polynomials.
(67)
and equation s is m
(1- ahJss)kl,- ah E J8jklj = hf,(yi) J=l
(68)
Ifs
Dividing by akl, and substituting using eq 66
and, rearranging
Occasionally first-order approximations for the derivative, eq 65, are insufficiently accurate, and we obtain a significant improvement using AYj AYa
+ w2k2j W l k l a + W2k% wlklj
(71)
Substituting in eq 64
Note that for real y,, & is real and when y, and dy,/dt are both positive or negative, i;, is positive. When ys and dy,/dt differ in sign, & is negative. Furthermore, for inherently stable ODES with real eigenvalues, when X, is negative, it approximates the largest eigenvalue, in m a g nitude, that contributes to the solution. This is illustrated in the results section. Furthermore, the stiff variable is not permitted to be a trace quantity, except when & is a trace quantity. This can occur is oscillating systems, where one variable passes through zero while another passes through a peak or a valley. Step-Size Adjustment. After estimating X,, error tolerances are checked and, when satisfied, the time-step counter i is advanced prior to step-size adjustment. Compared with the Michelsen algorithm (eq 30 and 31), we have developed a more efficient, but somewhat more complex, extrapolation algorithm well-suited for single-step integration formulas. Let y,. 2, yai-l,and yai be estimates of the second-order formu$ for the stiff variable and compute the coefficients for a quadratic polynomial through these points = P2t2 + P l t + Po
(76) This polynomial, illustrated in Figure 5, is used to extrapolate past ti. In addition, we let y’,j be the estimate of the imbedded, first-order formula, and assume the first-order estimate varies linearly with time. The slope of the straight line passing through ys,i-l and y;,i is Ya
Then, eq 68 and the sth row of eq 23
(77) To estimate performance of the first-order ODE past ti, a parallel line is passed through ysi with the equation Y i = P’lt
+ 8’0
(78) Equations 76 and 78 are combined to give estimates for relative local truncation error as a function of time ys - Y i r , = - - YS
1-h
x, =
Wlkls
ha
+ W2k%
(75)
- P 2 t 2 + (01- P’lN + Po - p’o
(79) Pzt2 + P l t + Po Then eq 79 is solved for t = t*i+lwith r, = e,: the relative error tolerance for the stiff variable. Equation 79 has up to four real roots, as illustrated in Figure 5. We use the smallest root that gives a positive step-size.
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981 261
Table 11. Ozone Decomposition at ti = 50. Adaptive Algorithm
a
E
no. of steps
Yl
% error
yz
% error
comput. time,a s
0.5 0.1 0.01 exact solution
30 61 98
0.000436 0.000431 0.000429 0.000428
1.9 0.7 0.2
0.014148 0.014003 0.013930 0.013872
2.0 0.9 0.4
0.6 1.2 1.8
UNIVAC 90/70computer, double precision arithmetic.
Normally, the candidate step-size, h* = t*i+l - ti, is a conservative estimate for a decaying variable because y s and y', straddle the true solution y s ( t ] ,whereas h* is an optimistic estimate for an increasing variable because y s and y', usually exceed y,(t). Therefore, we recommend h = ah* = a(t*i+l - ti) (80)
1. 2.
Set y { O > and s e l e c t h and Y.
-
Sat step-counter i Evaluate
0
giql
Caapute parameters a,b,w1,w2,W 5.
Compute
q+l,$,
Compute L1,L2
with
Caopute e + = w(k -k ) -2 -1
&O
a>l a 600 19 10.5 1.4 0.5
1.42 2.12 3.06 4.32 5.87
5.189 4.54
14.3
4.22
284
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981
Table VII. Eigenvalues for Linear System ( A , = -1, h , = -100, h 3 = -1014) time 1 x lO'l6 1 x 10-13 5 x 10-13 5x 2 x 10-9 1.5 X 2.5 x 10-4 0.014 0.022 0.036 0.060 0.072 0.089 0.119 0.183 5 9 a
stiff variable 3 3 3 3 3
2 2 2 2 2 2 2 2 2
-
-1014 -1014 -1 014 -2.2 x 1013 -2.6 X l o 8 consecutive steps -2.05 x 103 -100.02 -98.2 -96 -85.7
I
h 1.5676616 0.72859 X 101.5704603 0.27986 X 10-5 1.5707766 0.31632 X 1000t
1.5708087 0.32069 X 10-7 1.5708128 0.41630 X 1.5708128 0.83835
X
10-13
1 , 2, or 3 1, 2, or 3
Yz
no. of rejected steps
0.26936 X
0.85448
0
0.30127
X
0.85451
0
030961 X 10-4 0.35574 X 10-5 0.12558 X lo-" 0.18469 X 10-14 -094096 X
0.85452
0
0.85452
0
0.85452
0
0.85452
0
0.85452
0
X
0.85452
0
X
0.85416
1
X
0.84110
0
h*a
10-9
1 , 2, or 3
-0.997 -1.04 -1.01
First-order approximation, eq 70.
Table VIII. Locations of Points of Inflection of the Linear System with Complex Eigenvalues lOOOt stiff AA predicts true location variable 1.5708128 3.1416596 4.7126127 6.2835669 7.8544226 9.4252761 10.9961230
Table IX. Performance of the Step-Size Control Algorithm near the First Extremum (Peak in y 2 ) O
1.5707963 (n/2) 3.1415926 (n) 4.7123890 (3n/2) 6.2831853 (2n) 7.8539816 ( 5 s / 2 ) 9.4247780 (3n) 10.9955743 (7n/2)
1 2 1 2 1 2 1
algorithm, but many large systems can be decomposed into subsystems with similar rates of change, and elements of J need be evaluated only for the rapidly changing subsystems. We have demonstrated this approach for integration of the stiff ODES that represent the dynamics of a distillation tower (Prokopakis and Seider, 1980). Oscillating Systems. The linear system with complex eigenvalues is moderately stiff (SR = lo4). Yet, it involves a high-frequency oscillation that requires continuous adjustment of the stepsize through peaks, valleys, and points of inflection. Table VI11 locates the points of inflection in the region 0 < t I0.011 and Table IX shows performance of the step-size control algorithm over consecutive steps near the first peak in yz. In this region, y1 becomes less than but with y2 as the stiff variable, X, is less than lo4. Hence, y1 is selected as the stiff variable in Step 6 of AA. Careful examination of Table IX shows a sharp reduction of step-size through the point of inflection in yl. This is due to increasing rs as the curvature of the quadratic extrapolation reverses. Fortunately, the step-size control algorithm reduces h in a few steps and recovers as rapidly. It may be possible to avoid this reduction by examining the radius of curvature, although we have not yet undertaken such a study. Note that with a = 1.0 (& < 0) and a = 0.8 (& > 0), only six steps were rejected with 139 successful steps in the region 0 5 t 5 0.1. Whereas with a = 1.5 and 0.8,243 steps were rejected, with 287 successful steps. It should be possible for AA to adjust a to lower the number of steps rejected, although we have not explored this possibility. At t = 0.1 (after the exponential terms have decayed), with c = 0.1, high accuracy is achieved (ylwithin 0.006% and y2 within 0.005% of the true solutions) in 139 steps and 3.2 s of UNIVAC 90/70 computer time. Double precision arithmetic was used. Results for the Belousov Reactions are in Table X. For comparison, results are tabulated using Chan's IMB5 al-
1.5708128 0.i9245 X 10-l6 1.5708128 0.84567 X 10-1' 1.5749229 0.41101 X 10-5 1.6750097 0.10009 X 10-3 a
The values of
01
YI
10-3
10-14
-0.16636 10-13
-0.35107 10-2 -0.87892 10-l
are 1 . 0 (X, < 1 ) and 0.8 0,> 1 ) .
gorithm and the GEAR program. Since Chan and coworkers (1978) do not report their method of adjusting step-size, Eq 79 and 80 were used with (Y = 1. Chan and co-workers report 0.02% error for y1 at t = 300 with 288 steps and 6 = 0.01. We obtained 445 steps and 1.04% error, but with 6 = 0.005, we found 490 steps and 0.17% error. Clearly, AA is competitivewith IMB5 for integration of this limit cycle. The GEAR program performs especially well for 0 < t < 1.15, but as integration proceeds AA accumulates less error and for comparable accuracy at t = 300, computation times are comparable. Note that function evaluations consume s per step, and Jacobian evaluations 6 X lo4 s per step (several elements are constant). The remaining computations consume 145 X lO-'s per step for AA and 142 X s per step for the GEAR program. Figure 7a shows computed concentration profiles over a full cycle and Figure 7b expands the region 0 It I5 with the bold lines indicating the choice of stiff variable. In the region 0 It < 1.14, y1 is the stiff variable; followed by y2 for 1.14 < t < 1.166, y3 for 1.166 C t < 1.176, y2 for 1.176 < t < 3.486, y1 or y2 for 3.486 < t < 3.563, y1 for 3.563 < t < 3.846, and y2 from 3.846 < t < 5.0. Table XI shows the comparison of eigenvalues and pseudc-eigenvalues at various times during the integration. As expected, when the stiff variable decreases, & lies between the two largest IRe(Ajjl that contribute significantly to the solution; otherwise, when the stiff variable increases, & is positive. In this m e , the second-order approximation, eq 75, was necessary to give sufficiently accurate La. A reduction of h to 5 X loF4occurs in the region 1.146 It I1.166, where the variables change rapidly. Similarly, a reduction in h to 2 X lov4occurs in the region 3.907 I t I3.925, where y1 decreases by five orders of magnitude. Conclusions The following conclusions are reached. (1) The adaptive algorithm has been demonstrated to integrate accurately linear and nonlinear stiff systems with and without oscillations. (2) The rate of change of the stiff variable can be approximated by
dY,/dt = L Y E
(57)
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981 285
Table X. Performance of Adaptive Algorithm for Belousov Reactions ~______
E
t = 1.15 no. funcof Jacobian tion steps evals. evals.
%error
t = 4.00 no. funcof Jacobian tion steps evals. evals. % error Adaptive Algorithm 245 245 516 1.02 361 361 744 0.04 665 665 1336 0.0
0.01 0.005 0.001
115 160 238
115 160 238
242 324 476
1.33 0.06 0.05
0.01 0.005
138 144
138 144
288 292
1.28 0.36
367 395
97 139 199 289
14 14 18 23
169 215 278 365
0.3 0.05 0.007 0.001
280 391 520 716
10-5 10- 7
t = 300.00 funcno. of Jacobian tion steps evals. evals.
time,a %error
5
342 434 812
342 434 812
740 912 1636
1.20 0.06 0.002
9.2 11.6 21.8
Chan’s IMB5 367 756 395 810
1.01 0.08
445 490
445 490
936 1016
1.04 0.17
13.0 14.5
GEAR 495 602 751 938
1.47 1.11 0.08 0.07
391 556 752 1035
59 64 77 98
684 923 1219 1548
9.01 2.38 0.14 0.02
5.5 7.6 10.2 13.6
41 48 49 56
UNIVAC 90/70 computer, double precision arithmetic. per step.
(Y
= 1in eq 80. Requires additional matrix multiplications
Table XI. Eigenvalues for Belousov Reactions
t 0.1969 1.057 1.137 1.145 1.152 1.164 1.168 1.243 2.990 3.864 3.916 3.921 3.994 14.24 45.92 64.41 115.7 219.6 286.1 292.8 298.9 300.0
h3
A1
-84.082 -60.565 -7.949 -16.339 -25.381 -51.985 -64.401 -152.20 -104.70 -36.480 -612.06 -771.13 -3018.6 -1.33 X 10’ -7.76 x 104 -48245 -12807 -868.36 -152.62 -126.25 -104.03 -99.744
-0.0933 * 0.1172i -0.1020 * 0.96843 -112.65 -1.6311 -178.20 -0.9162 -261.43 -0.6524 -523.04 -0.40 10 -646.66 -0.3546 -1524.8 -0.2428 -1100.5 -0.2772 -509.61 -0.4521 -0.0934 f 0.0759i -0.0287 -0.1582 -0.0259 -0.1610 -0.0259 -0.1610 -0.0259 -0.1610 -0.0259 -0.1610 -0.0259 -0.1610 -0.0260 -0.1610 -0.0344 -0.1 525 -0.0429 -0.1440 -0.0825 -0.1043 - 0.0934 f 0.0310i
stiff variable Y1 Y1 Y1
Yz Yz
Yz Y3
Yz
y2 Yz Y1 Y1 Y1 Y3
Y3 Y3 Y2
Yz Y1 Y1
Y3 Y3
AsU
-2.954 32.257 71.108 -104.211 -190.039 -551.617 44.555 1.973 -1.626 41.794 -47.437 -644.1 77 -0.0187 -0.161 -0.161 -0.161 -0.0259 -0.0259 -0.3547 -0.0699 0.1965 0.2611
stiffness ratio 901 594 69 195 401 1304 1824 6280 3970 1127 6553 26869 116548 5.14 x 107 3 x lo6 1.86 X l o 6 494479 33398 4437 2943 1261 1068
Second-order approximation, eq 75.
Using kl,8 and k2,,over the current time-step, the pseudo-eigenvalue, &, can be estimated and used to adapt the parameters of a second-order, semi-implicit Runge-Kutta method to give improved accuracy for integration of the stiff system over the next time-step. (3) The three most recent values of the stiff variable, ys,i,ysi-l, and ys,i-2 (high-order estimates) and y;,i (firstorder estimate) provide an excellent basis for step-size adjustment. The algorithm has been demonstrated to adjust step-size by ten orders of magnitude in five steps. It has successfully traced System 4 which oscillates with high frequency. Some inefficiency has been observed in passing through points of inflection, where for System 4 the step-size underwent a sharp reduction in a few steps followed by a rapid increase. (4) When adjusting step-size, the parameter a should be set to avoid large numbers of rejected steps. With a = 1.5 (X, < 0) and 0.8 (X, > 0), few rejected steps occur for most systems. However, for System 4, 243 steps were rejected, whereas with a = 1.0 and 0.8, only six steps were rejected. (5) For nonoscillatingsystems, when dy,/dt and & differ in sign, X, gives a fine approximation to the largest ei-
genvalue, in magnitude, that contributes to the solution. (6) In limited tests, the adaptive algorithm is competitive with the two most efficient semi-implicit Runge-Kutta algorithms reviewed by Chan and co-workers (1978): IMB4-Michelsen’s third- and first-order imbedded pair, and IMB5-Chan’s fourth- and second-order imbedded pair. In some cases, the adaptive algorithm is more efficient. (7) In limited tests,the adaptive algorithm is competitive with the GEAR program, Revision 3 (Hindmarch, 1974). GEAR does not evaluate the Jacobian each step, but takes more steps than the adaptive algorithm for the stiff linear system 3 and for the Belousov reaction system 5. (8) The adaptive algorithm has been implemented as a prototype program in approximately 200 FORTRAN statements. It is compact and well-suited for process control computers. Acknowledgment The helpful comments of C. W. White, 111, are appreciated. Computing time was provided by the School of Engineering and Applied Science at the University of Pennsylvania and is acknowledged.
Ind. Eng. Chem. Fundam., Vol. 20, No. 3, 1981
266
Nomenclature a, al, u2, ai, = Parameters in semi-implicit Runge-Kutta methods b, bl, b2, b3 = parameters in semi-implicit Runge-Kutta methods b31,bS2= parameters in semi-implicit RungeKutta methods bdl, bd2, bM = parameters in semi-implicit Runge-Kutta methods c, c1 = parameters in semi-implicit Runge-Kutta methods A, B, C,D, E = constants in eq 61 C1 = (Ky1(01)/(X1- 1 2 ) in eq 6 C2 = y2M.- ( K Y I K $ I ( ~- X2) in eq 6 dj = jth eigenvector in eq 3 el?_= local truncation error for variable j f - function h = step-size I = identity matrix J = Jacobian matrix kj = increment functions in semi-implicit Runge-Kutta methods m = number of ODES q = see eq 30 and 31 r = order accuracy of semi-implicit Runge-Kutta methods rs = relative truncation error, eq 79 t = independent variable, time t*i+l = suggested final time for ith step u = number of terms wj = weighting factors in semi-implicit RungeKutta methods y = dependent variable ys = stiff variable-has largest fractional rate of change boldface italics = vector of dimension m boldface roman = matrix of dimension m X m Greek Symbols a = multiplying factor for step-size aj, Pj
= parameters in multistep methods; eq 32 and 34
&, PI, b2,Bo,B1= parameters in eq 79 for step-size adjustment y =
characteristic root of difference equation
= absolute error tolerance for variable j trj = relative error tolerance for variable j E,
= parameter in ozone decomposition system = parameter in ozone decomposition system & = eigenvalue (= XR + iXI) X, = pseudo-eigenvalue-defined in eq 57 t = value of dependent variable ul,u2, u3 = parameters of linear system with real eigenvalues Subscripts 0 = initial value at t = 0 a = variable at -hX = max = maximum E
K
min = minimum s = stiff variable Superscripts - = pseudo* = suggested ’ = first-order estimate
Literature Cited Alken, R. C.; Lapidus, L. AICE J . 1974, 20, 388. Ballard, D.; Broelbw, C. “Dynamic Simulation of MutDieWtbn Columns”; presented at the 71st Annual Meeting of AICM, Mlami, Nov
1978. Bowen, J. R.; Acrlvos. A.; Oppenheim, A. K. Chem. €ng. Scl. 1963, 78,
177. Byrne, G. D. prhrate communication, 1981. Byme, 0.D.; Hlndmarsh, A. C. “Appilcatkns of EPISODE. An Experbnentai Package for the Integration of Systems of ODE”, in ”Numerical Methods for Dlfferentlel Systems”;Lapldus, L., Schlesser, W. E., Eds.; Academic Press: New Ywk, 1976,p 147. Calllaud, J. B.; Padmanabhen, L. Chem. €ng. J. 1971, 2, 227. Calahan, D. A. Proc. IEEEl966, 56, 744. Cash, J. R. ACM J. 1976, 23, 455. Chan, Y. N. I.; Bknbaum, I.; Lapidus. L. Ind. Eng. Chem. Fund“ 1976,
77, 133. Ehle, B. L. ”A Comparison of Numerical Methods for Sohrlng Certaln Swf OM”;Vlctorle Unhrersity, Report No. 70, 1974. Enrlght, W. H.; Hull, T. E. “Comparing Nwnerlcal Methods for the solution of Stiff Systems of OMSArising in Chemistry”, in “Nwnerlcel Methods for Differential Systems”; Lapidus, L., Schlesser, W. E., Eds.; Academic Press: New York, 1978,p 45. Fleld, R. J.; Noyes, R. M. J. Chem. phys. 1974, 80, 1877. Gear, C. W. “Numerical Initlei Value Problems in Ordlnary Differentlei Equations”; RentlceHell: Englewood Cliffs, NJ, 1971. Hlndmarsh, A. C. “GEAR: ordinery Differential Equatkn System Solver”, Computer Doclgnentetion, UCID-30001,Rev. 3,Lawrence Livermore Lab oratory, Dec 1974. Hindmarsh, A. C.; Byme, G. D. “EPISOM: An Effective Package for the Integration of System of OrBnary Differential Equations”, Computer Documentation, UCID-30112, Rev. 1, Lawrence Llvermore Laboratory, April
1977. Lapidus, L.; Seinfeid, J. H. “Numerical Solutbn of Ordinary Dlfferential Equations”; Academic Press: New York, 1971. Michelsen, M. L. AICM J. 1976, 22, 594. Prdtopakis, 0.J.; Sekler, W. D. “Dynamic Simulation of DistillationTowers”; presented at 73rd Annual Meeting of AICM, Chicago, Nov 1980. Robertson, H. H. “Solution of a Set of Reaction Rate Equations”, in ”Numerical Analysis”; Walsh, J., Ed.; Thompson Book Co.: Washington, 1967,p 178. Rosenbrock, H. H. Comput. J. 1863, 5, 329. Seinfeld, J. H.;Lapidus, L.; Hwang, M. Ind. €ng. Chem. Fundem. 1970, 9 ,
268. Shampine, L. F.; Davenport, S. M.; Watts, H. A. "Comparison of Some Codes for the Initlai Value Problem for OM”, “Numerical Sokrtkns of Boundary Value Problem for WE”; Azlz, A. K., Ed.; Academic Press: New York, 1975,p 349. Viiladsen, J.; Mlcheisen, M. L. “Sokrtkn of Dlfferential Equatbn Modeis by Polynomlal Approxlrnatbn”; PrentlcsHall; Englewood CimS, NJ, 1978.
Receiued for reuiew July 18, 1980 Accepted March 31, 1981
