- Email: [email protected]
Solving STIFF systems by Taylor series
251
Solving STIPF systems
by Taylor
Series
Y. F. C-hang 976 U
Foothill BLvd.. Claremont. CA. 91711
Abstract The analysts ot long Taylor series for the solutions of stiff prohlems gives new insight into their true character and leads to new solution algorithms. Some ODE’s are difficult to solve becau se of singularities close to the real axis. This is not true for stift problems. There is a linear ODE hidden within the stiff system of ODE’s. When the hioden linear ODE is of order one, there is an exponential component in the solution. Most stiff problems arc of this type. This exponential forces the integration stepsize to be very small relative to the domarn of the stiff problems. This exponential component is readily made apparent by analysis of the Taylor series for the solutions. We analyze in detail a stiff problem from chemistry, which is very difficult to solve. The portions of the equations giving rise to the stiffness in this problem are identified. In this problem. there are more than one hidden linear ODE. Artificial stiff problems with up to stiff-order four have been created by combining linear ODE’s with a non-linear problem with known solutions. The Taylor series algorithm for solving stiff problems is an integral part of the ATOMFT system for the automatic solution of ODE’s. Introduction:-
Taylor
Series Methods
The use of Taylor series to solve ODE’s began with Sir Issac Newton, who of solved ODE’s with 4-term Taylor series. In modern times, the development Taylor series methods for ODE’s was independently pursued by Gibbons(14), Moore(21), Leavitt(“), Barton, et al(l), Chang, et al(3*8), Rail(“), and Halin(15). They were all able to demonstrate that Taylor series methods do make a practical algorithm for the solution of ODE’s. This author’s Taylor was shown to be faster than the usual series package called ATOMCC(*) methods when solving the problems in the Toronto non-stiff ODE test set at high precisions. There is at leact one book (**I that advocates the Taylor series methods as general purpose algorithms for ODE’s. There are at least four misconceptions concerning the efficiencv. stability, and general. viability of the Taylor series methods as applied to ODZ’s. They ale:- (1) difficulty to evaluate the high order derivatives, (2) sequential differentiation by the use of recursion relations is grossly unstable, (3) Henrrci’s ‘proof’(16), and (4) a ‘proof’ that the high order derivatives may be unbounded. The facts speak!: against these misconceptions. Analyses of the high order Taylor series terms in ODE problems yield the precise location and the precise @ELsevier Science Pubhshing Co., Inc., 1989 655 Avenue of the Ameruxs, New York, NY
10010
0096-3003/99/$03.50
252
order of the singularities in the solution functions. Such precision would be impossible without numerical stability. We shall answer these misconceptions in the following paragraphs. ( 1)Many numerical analysis texts echoes the incorrect conventional wisdom that Taylor series methods are impractical due to difficulties in computing the Longs Tay&r2;e;iy are sequentially evaluated using high order terms. with accuracy and efficiency. It is not differentiationarithmetics - * a svmbolic process. and it does not use finite differences. There is a recent hibiiography(ll) containing 90 references to differentiation arithmetics and its applications. The recursion relations used in differentiation arithmetics are derived from the OPE’s themselves by applications of the chain rule. There is no difficulty in evaluating high order Taylor series terms. l
(2)There are those who hold to the view that the sequential evaluation of Taj-lor series terms is grossly unstable and leads to unusable series terms. Although it is true that all sequential processes are inherently unstable, this is not of any major concern in the evaluation of Taylor series terms, because the magnitude of the instability is usually very small. Except for one class of problems. w-e have found ail recursion relations to yield very good results. A Taylor series of 10.000 terms have been calculated with high precision for the eIliptic-integral function. The one class of problems, where the instability does cause some concern, are called ‘last-pole’ problemsC5). This error condition occurs when the solution in these problems approaches an aysmptotic entire function, which does not have any singularities in the finite complex plane. The difficulty arises as the solution passes beyond the ‘last-pole’ in the non-linear domain of the problem into the linear domain. There are subtractive errors that arise from the ODE’s becoming almost linear with a vanishingly small (and so the subtractive error) non-linear component. (3)It was unfortur.ate for Henrici(16) to have offered a proof that analytic continuation obtained by rearranging a truncated Taylor series is unstable. Although his proof is correct, it is totally irrelevant. None of the Taylor series methods developed for ODE’s performs analytic continuation by rearranging the terms of a truncated series. They ail generate a new Taylor series for every function at each new integration step. They all have successfully solved ODE’s through hundreds of integration steps extending far beyond the domain of Henrici’s concern. (4)There is another ‘proof’ that is doing the rounds of numerical analysis texts. It deals with the fact that high order derivatives increase without bound when there are singularities close by. This authorC4) has shown this proof to be correct but also totally irrelevant. The derivatives do increase without bound; however. the Tqlo: series terms are bounded. As a matter of fact, the Taylor series terms typically describe a damped sine wave. Figure 1 is a graph of the first 100 terms of F(n) versus n at x=0 for the solution off” + f’f + f = 0, withf(O)=-2.81 andf’(O)=-0.99. There are simple poles at 2=2.2 f i0.69; the stepsize used (h=2.3) is slightly less than the radius of convergence. The angle of the poles measured from the real axis at x=0 (17.5”). is exactI) equal to the phase angle difference between the Taylor series terms.
Figure 1. Typical Variation
of Taylor Series Terms from ODE’s
The Taylor series methods implemented by each of the Barton, Chang. and Halin teams were compilers that accepted mathematical statements of the differential equations to be solved. The compilers produced an object program that is executed to solve a specific problem. Thus. problems are solved ‘tiatomatically’: no formal programming is necessary. This represents a significant advantage over the usual methods. The original compiler(3) written by this author was computer specific. it Further was superceeded by a portable program called ATOMCC(4*s’. extensions to include the solution of stiff problems, arbitrary functions definable by the user, and the soiution of DAE’s have led to the present ATOMFTc6) program. The Taylor series method 4as evolved into a general purpose numerical technique capable of solving all manner of problems. The ATOMFT packagec6) applies the detailed analyses(3*7) of the Taylor series terms to advance the understanding and treatment of the truncation Tne Taylor series analyses yield the location errors in numerical integration. and order of the singularities in the solution functions. Thus, ATOMFT is able to control local errors with sufficient accuracy to achieve proportionally between the local trunation and the global errors. This also represents a significant advantage over the usual numerical integration techniques. Introduction:-
Taylor Series Solution of Stiff Problems
Stiff s stems of ODE’s were first discussed by Curtiss et alc9’ in 1952. as the source of the difficulty Dahlquist(“) had identified numerical instabilit conventional Enright, et al (12.13 examined in solving these problems. numerical methods for stiff problems. The Taylor series method to be discussed is different from conventional methods. First, we use 30-term Taylor series; thus, it is a 29th order method. Second, we analyze the series to obtain the radius of convergence and other characteristic constants such as exponential coefficients. Third, the integration stepsizes used are much larger than in other methods and they are infinitely variable.
254
In 1967, while developing the analysis of Tayior series, we were given a problem by Barnes, et al CL). This problem was stiff. It did not obey any of the Recent studies has revealed simple rules or formulae being then developed. recursion relationships identifiable with stiff problems. The numerical instability inherenr in the recursively generated Taylor series terms does I.ot show up in any casual inspection of the Taylor terms. In contrast to the last-pole problems mentioned above, most if not all non-linear problems have Taylor series terms that are reliably reproduced. Rather than the propagation of large subtractive errors as in last-pole problems, the numerical instability in stiff problems is quite difficult to detect. We need to ha-fe a clear understanding of this instability in order to accurately and reliably evaluate an exponential coefficient A. The singificance of this exponential coefficient is explained in our model stiff problem in the next section. In our solution algorithm, once the coefficient A is found to be constant, we truncate the Taylor series to remove ALL of the terms that are numerically unstable. Therefore, we have the best of both worlds. The numerical instability is analyzed to find a constant coefficient A. Then, that value of x is used to take giant integration steps, while all those Taylor terms used to evaluate that x are removed. A Model Stiff Problem
When long Taylor series for every function in a problem is available for analysis, simple relationships between adjacent Taylor terms become readily apparent. Assume, for the moment, that a stiff Frcblem contains a first order linear ODE within a system of non-linear ODE%. This assumption will be justified by the analysis of data from actual stiff problems. Therefore, we have
whose solution is the exponential function, s(t)=exp(-At). the total solution for the stiff problem is of the form
f(t) wheref(t)
=
S(t)
+
g(t)
We also assume that
,
is the total solution, and g(t) arises from the non-linear
[21 stiff ODE%.
The adjacent Taylor terms of the S-series are related by:-
*
(
n-
l)e+=-x,
where n is the order of the terms, and h is the stepsize. When the magnitude of the exponential coefficient, x, is large, the magnitude of S(n) increases very rapidly with n. The magnitude of S(n) reaches a maximum at some very large value of n, and then decreases for all higher n’s, The peak of this ‘hump’ for the
usual stiff problem occurs beyond beyond our sphere of interest.
hundreds
or even thousands
of terms,
far
The magnitude of the Taylor term, G(n), is determined by the positions of the singularities in the solution in the complex plane of the independent variable t. The radius of convergence is defined to be the distance from the point of expansion to the nearest singularity. If the stepsize is equal in magnitude to or less than) the radius of convergence, the Taylor term, G(n), is bounded(5* s 1. Typically, the Taylor term, G(n), describes a damped sine wave (see Figure 1). This is true for solutions of -all non-linear ODE systems. Since the Taylor series for g(t) is bounded. and the function s(t) has terms that increase very rapidly with n, the high order F(n) are dominated by S(n). Therefore, if x is large, we cau expect the same relationship between adjacent Taylor terms of f(t) as for s(t) when n is large, say around IO.
P- 1)y+--a (
n-
When A is large, the ratio F(n)/F(n-1) is large except when the stepsize is very small. This is the source of the stiffness in stiff problems, based on the assumption stated in Eqs.[ I] and [2]. When solving a non-stiff non-linear ODE system, the stepsize is restricted to be less than the radius of convergence. If there are singularities close to the real axis, the solution must proceed with small stepsizes. Any attempt to take large steps will result in solutions that have very large errors. In contrast, the solution of a stiff system of ODE’s has a stepsize forced to be small by an exponential component. In this case, taking giant steps should not result in appreciable errors if we account properly tnr the singularities in the non-linear g(t) function. This is our model of stiff prol!tms. We will attempt to fit this model to some actual stiff problems. The stiff-order of a stiff problem is defined to correspond to the order of the hidden linear ODE. Almost all of the non-linear problems in the EnrightHull set(13* 14) with and without the modifications recommended by Shampine(17) have stiff-order one. We will examine stiff problems with stifforder higher than one after completing the analysis to find the exponential coefficient x.
Detection
of Exponentials,
Stiff-Order
One
In this section, the analysis of Taylor series problems in the search for a reliable evaluation of the The numerical instability discussed by Dahlquist(“) We have found it necessary to analyze this instability evaluate the exponential coefficient x.
is applied to many stiff exponential coefficient X. is very much in evidence. in full in order to reliably
Barnes(l) problem was concerned with the production, inside a run-away nuclear reactor, of radioactive methyl-iodide (CH3I) from methane (CH4), free
iodine (12). hydrogen iodide (HI). the iodine ion (I). and the methyl ion (CH3). The iodide reaction equations for the nuclear accident problem of Barnes are:n I2 _-.
dt
- - kj 12 # + k6 1' H .
d CH4
- dt
d CH3I -__--dt
_ k3 CH3 HI - k, CH4 I ,
141
= _ k: CH31 I + kZ CH3 I2 .
M = I + CH4 + CH3 + HI + 12 + CH31 .
with the initial concentrations
given by
2 12(O) = 2 12 + HI + I + CH31 . .G CH4(0) = 4 CH4 + 3 CH3 + HI + 3 CH3I . and CH4(0) = CH4 + CH3 + CH31 .
The initial concentrations
are:- CH4(0)=4x10‘8.
The reaction rate constants. k’s, depend nuclear reactor. They are given in Table 1. TABLE Temperature
%% 600K
kl
8.31E+6 6_30E+5 8_7OE+3
1. The Reaction kz
4.67E+9 3.89E+9 2.85E+9
12(O)= 10‘lo, and CH31(0)=0. on the temperature
Rate Constants. ks
ks l_OOE+9 7.43E+8 4.57&+8
inside the
2.46E+4 3.55E+2 0.317
kg
5.51E+4 1.43E+9 ;.9887v2 1.92E+9 2.84E+9 .
The numerical solution of this problem at the temperature of IOOOKis quite easy. The problem is hardly stiff at all. and it is easily solved by most numerical methods. At the temperature of 80OK, this probiem is stiff but not difficult to solve. At the temperature of 6OOK. this problem is extremely stiff. The chemical reaction times to reach steady state inside the nuclear reactor are:- 10,000 seconds for LOOOK,300.000 seconds for 80OK, and 60,000,OOO seconds for 600K. These are chemical reaction times and NOT computer calculation times. Under the assumed model stiff problem of Eqs.[ I] and [2]. the ratio given in Eq.[3] should yield a coefficient A that is constant. We analyzed the Taylor series for 12 in the three temperatures for the Barnes problem. Listed in TABLE 2 are values of x for all Taylor orders from n=7 to n=30 at a fixed time. These ratios are calculated with double-precision; therefore, they should be accurate to at least 12-digits. For our proposed model of stiff problems to be valid, the
257
numbers listed in TABLE 2 should be constant the Taylor order, n. TABLE 2. The Exponential n
lOOOK :
BOOK
.389177 “4:29:3: .389214 .389250
:%26d: .389286 .485546 .487678 .489805 491926
: %Ei
:"499% .498290 .500440 .5026?9
.I%9500 :388%%
: 5500% .510336
-389616 .389623 .389601
.513975 .519020
.526669 .539000 .559320 .592065 :sM;
: %k% .389465
.389595
: %53% .388892 .387993 : %4% .374558 .357958
in each column independent
Coefficients, 600K
of
x. Krogh
1000.03 999.99 999.99 999.99 Ext 999:92 999.84 999.69 999.38 E% ;;;:cI'j 979:a7 958.93 914.36 812.71 539.24 -5E ii 2654.85 2244.88 2106.70
In fact, the x’s calculated are not ctinstant to the degree expected, casting doubts on the validity of the model. Further study is needed to establish the truth of this model. Also, if this model is correct, the stiffness is not determined by the magnitude of the exponential coefficient, x. Most of the coefficients in TABLE 2 are smaller than 0.5. The 6OOK problem is the most stiff among the three, and yet it has the smallest coefficient. Is the stiffness of a problem indicated by the constancy of the coefficient. 3 This is shown to be untrue by the results [in the last column of TABLE 21 for the very stiff Krogh problem(18) (also E4(12)). The calculated coefficients for Krogh’s problem are certainly not constant. Dahlquist(l’) had identified numerical instability as the source of the it behooves us to suppose that difficulty in solving these problems. Therefore, the non-constaucy of x shown in TABLE 2 is due to an accumulation of errors. For the moment, let this error propagation be linearly dependent on the order of the terms. We will fit the data by a straight line x,
=
x,
+
an,
161
where X, is the ratio listed in TABLE 2, x0 is this ratio extrapolation to the zeroth term (the vertical intercept), n is the order of the Taylor term, and o is the slope. Results of applying this fit to the Barnes problem and many other stiff problems indicate that the numerical instability is more complex than a simple linear propagation of errors.
258
We next fit the non-constant x in TABLE 2 to a propagation of errors that is plus an also dependent on the power of n. This fit is a linear correction exponential correction; &I
=
.b
+
an
+ ar”
i7]
.
where 6 is a constant, and $’ is the exponential component With this approximaticn, all the stiff problems examined constant with $-digit precision.
of the correction. yield A’\ that are
The results for the Barnes problem at 80OK are listed in TABLE 3. The first three terms cannot be fitted because the linear exponential function. s(t) has yet to be dominant. The next six terms have a linear behavior. All the remaining terms required the full approximation. These results are typical of all the stiff problems examined. TABLE 3. Calculated n 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
raw X
error
-3.8917D-03 4.3D-02 -3.8921D-01 9_6D-05 -3.89250-01 9.2D-05 -3.89280-01 1.8D-07 -3.8932&01 -3.89351)-01 -3.89398-01 -3.8942D-01 -3_8946D-01 -3.895OD-01 -3.89531)-03j-g:;; -3.8956D-01 7:aD-oa -3.8959D-01 8.4D-08 -3.a961D-01 -3.89621)-01P*s% -3.8960D-01 2:5D107 -3.8952D-01 8.2D-07 -3.8932D-01 3.38-06
Exponential fit X
Coefficients.
Q
-3.88998D-01 -3.60D-05 -3_88998D-01 -3.59D-05 -3_8899aD-01 -3.59D-05 -3.889991)-01-3.598-05 -3.89000D-01 -3.57D-05 -3_89003D-01 -3.55D-05 -3.88998D-01 -3,60D-05 -3_88998D-01 -3.6OD-05 -3.88998D-01 -3_6OD-05 -3.88998D-01 -3.60"-05 -3.88998D-01 -3.608-05 -3.88998D-01 -3_60D-05 -3.88998D-01 ->.60D-05 -3.8899aD-01 -3.59D-05 -3.89000D-01 -3.59D-05
80OK.
B
g:;:
-7
z%z:
17D-10 z:oD+oo 15D-10 2,oD+oo z.oD+Oo gr,l:: 2.oD+oo 13D-10 12D-10 EE% llD-10 2:oD+oo
The important conclusion of this study is not our ability to the fit the data with a complex error analysis; it is the clear existence of a constant A. In the fourth coIumn, the value of A &es by 5 in the sixth digit. Therefore, we conclude that the model for stiff problems, Eqs.[l] and [2J, is valid. A problem is stiff if upon application of Eqs.[6] and ]7] a constant A is found. Incidentally, the accuracy of x does not materially effect the accuracy of the ODE solution. As is explained later in the section on the stiff algorithm, the exponentia1 function, exp(-at), is first subtracted from the solution and then added back at a later stage. It bears repeating that, in our solution algorithm, once the coefficient x is found to be constant, we truncate the Taylor series to remove ALL of the terms that are contaminated by numerical instability
259
Synthesis
of a Stiff Problem
In this section, we further demonstrate the validity of our model for stiff problems by synthesizing a stiff problem. The simplest model of a stiff problem consists of a system of non-linear ODE’s for the function g(t) plus a linear ODE, Eq.[ 1). Hence, it should be possible to create a stiff problem artificiallv. Until we can create a stiff problem out of known components, we cannot claim to understand stiff problems. This task is undertaken equations can be as follows:-
c&I
dt
-
p
s!.E ’
dt
=
in small steps.
6 g”
+
t
+
s ,
First, a possible
stiff system of
and
with g(0) = s(0) = 1, and p(O) =0 on [0, 11. The solution to the first two equations, Eq.[8] without the s at the end, is the first Painleve transcendent. To this well known problem, we add the function, s, and the associated linear ODE. The first Painleve transcendent, g, is a monotonicaly increasing function with a simple pole on the real axis at ? = 1.2. The addition of the linear function, s, to the Painleve transcendent does not As in any non-stiff solution, s(t), rapidly make this system a stiff problem. decays to some negligible magnitude for a large x. Therefore, beyond a very small step in time, the contribution of s(t) to the solution g(t) becomes negligible, and the solution proceeds ignoring the exponential. Therefore, this simple combination of ODE’s does not produce a stiff problem. In a stiff problem, as defined by our model, insignificance. It must somehow persist in disturbing the in spite of its rapid decay. Therefore, there must be a the non-linear ODE’s for g(t) and the linear ODE magnitude for s(t) can remain large. To accomplish replace the simple linear ODE by ds dt
= - A s
+
c g(t)
,
s(t) cannot decay to solution for large times cross-coupling between for s(t) such that the this cross-coupling, we
191
where c is the cross-coupling constact. This system of ODE’s is indeed stiff for non-vanishing values of c. We have succeeded in creating a stiff problem from known components. The stiffness of this problem is controlled by the magnitude of the coefficient, x, and the amount of the cross-coupling. After studying the results from this synthesis model, we conclude that the magnitude of x does not by itself determine the stiffness. Some very stiff problems have x’s that are quite small. The magnitude of the cross-couplmg constant isvery important in the stiffness of a problem.
260
Analysis of Stiff Problems In this sectioo, we take the final step towards a clear understanding of stiff problems. Since the model for stiff problems contains a linear ODE which has a cross-coupling component from the non-linear ODE’s, it should be possible to ideniify and isolate the source of stiffness in problems. It is important in this identification to remove all non-linear portions of each ODE to see whether it satisfies the form of Eq.191. The coefficient A is allowed to be a slowly-varying function of time. The final test, to identify the source of the stiffness in the problem, is to compare the value of the slowly varying function against the value of a calculated from the Taylor series. Of al! the stiff problems known to us, the most difficult ones are the Barnes methyl-iodide problem and CHEMP in Enright-Hu11(13). We now identify the sources of the stiffrzss in CHEMO. This problem has two distinct stiff phases, separated by phases of uormal behavior. The ODE’s for CHEMP are:-
ax
dt
= 77.27(y - xy + x - 8.3758-6 2)
xy + z)/77.27 ,
Q=(_y_
df ZL
,
and
= 0.161(x - 2) .
with initialvalues
x(O)-4,y(O)=
1.1, r(0)=4.
on [0,300].
To identify the source(s) of a problem’s stiffness, we look for occurrences of linear ODE’s. It is important to allow the exponential coefficient, A, to be a function of t. Only after having solved the problem, and having evaluated that function. can we determine if that function is nearly constant. The procedure is to force the ODE to fit the following form, collecting all the relevant portions into the X(t) and all the irrelevant portions into the [whatever]. ds zG=
- A(t) s
+
[whatever]
With three ODE’s in CHEM9, there are three potential model. They are:(1)
dx ;ii:=
(2)
d-1!= dt
-
(3)
dr z=
- (0.161) z
After
candidate,
- 77.27 (y - 1) X
$?$.
y
+ f
+
occurrences
of the above
[whatever] ,
[whatever] .
and
[whatever]
having solved CHEMY, we can eliminate ease(3) as a potential because at no time in the solution is there a x=0.161. That leaves
261
both ca\e( I) and cu\e(2) a\ possihle contenders for the source of the stiffness. Case( I ) has x(t) = 77.27 (v- I ). Therefore. .v(t) must he slowly varying for case( 1) to he the source. Case(Z) ha> x(t)=(x+ 1)/77.27. Therefore, x(t) must he slowly varytng for CHEMO to he stiff. So. the search is for phases of the solution when either u(t) or l:(t) I\ varying slowly. The complete
of CHEMV i\ shown in Figure 2. There is a sharp The vartation of y(t) is somewhat more gradual. The two functions are shown together in Figure 2: however, they differ in magnitude hy a factor of 70. The spike in x(t) is 70 times larger than the peak of r(t). An enlarged view of the solution near the origin is shown in Figure 3. spike
in .t(t)
solutton
near the origin.
linear -- vertical scale *. 9 .
.
oocl
.
.
-
x(t)
-
Y(t)
-. -. -. -... --.__ . ...**. 0
100
0
Figure 2. Solution of CHEMY. for x and y.
There are four distinct phases in the solution of this problem. The four phases are in the times [0, 1.256). (1.256.3 g), (3.8.3.98), and (3.98.300). 1. In the first phase IO. 1.256), our Taylor series solution shows a conjugate pair of singularittes at z=-0.35 + 10.30 and another pair at z= 1. I5 + iO.03. The second singularity pair causes x(t) to increase exponentially from x(0)=4.0 to x( 1.256)= 1.1784~10~. Here, CHEM9 is not a stiff problem. Whenever there is a rapidly varying function in a problem, that problem cannot he stiff. This is because a rapidly varying function is a clear indication that the solution is under the influence of some singularities. Then, the Taylor series of the solution cannot yield a constant value for X. As shown in Figure 1. the Taylor terms describe a damped sine wave when there are singularities close by.
262
In the second phase (1.256.3.8) x(t) is slowly varying. Therefore, ease(2) is the source of the stiffness. The calculated A agrees with the predicted value of (x+ 1)/77.27. The magnitude of x(t) decreases linearly from x(1.256)=1.1784x105with~=1.500tox(3.8)=4.968x104withx=692. 2. In the third phase (3.8, 3.98), x(t) decreases very rapidly from x(3.8)=4.968x10* to x(3.98)= 1.04. CHEM9 is again not a stiff problem. 3. In the fourth phase (3.98.3001, y(t) is slowly varying. Therefore. case( 1) is the source of the stiffness. The calculated x agrees with the expected value of 77.27 (y-l). The magnitude of y(t) increases slowly from y(3.98)=29 with x=2.2x103 to a maximum of y(17.1)=1.757xld with x = 1.36x105, and then decreases slowly to ~(300) = 1.22 with x = 17.
t
linear vertical scale stiff Case(P) stiff Case(l)
000
-
x(t)
. . . - Y(t)
5
6
7
a
9
10
Figure 3. Solution of CHEM9, an Enlarged View. The regions of stiff solutions are marked in Figure 3. It is clear in this graph that the spike in x(t) has a definite structure. The first stiff region does not begin until after x(t) has reached its peak and is decreasing on a sloped plateau.
263
The first stiff region ends when x(t) drops too rapidly !o sustain a ‘constant” x for Case(Z). The second stiff region begins once y(t) has becomed large enough to make the x for Case( 1) a dominant factor. The stiff algorithm is ahsolutly necessary in the solution of CHEMY. In the first stiff region. the integration stepsize!, are about ten times larger than a normal Taylor series stepsize. In the second stiff region. the integration stepsizes are at least a thousand times larger than the normal stepsizes. We now attempt to identify the stiffness in Barnes problem. Refering to the ODE’s in Eq.[4], if 12 is to be the linear function. then (M) must be slowly varying, but it is not. If CH4 is to be the linear function. then (HI) must be slowly varying. but it is not. If CH3I is to be the linear function. then (I) must be slowly varying, but it is not. Therefore, we must look elsewhere. After considerable following substitutions.
search, we found CH31 to be the linear function From Eq.[SJ, we identify CH3(t) to be
by the
CH3(t) = CH4(0) - CH4(t) - CH3I(t) .
Substituting this result into the term (k2 CH3 t2) in Eq.[7] and isolating the relevant and irrelevant parts, we find the linear first-order ODE that is the stiffness of this problem. d CH31 -.X dt
- (k,I2)CH31 +
[whatever],
where A(t)=k212. During the extensive SPYC~. we needed the availability of the solution, where it was determined that I2(t) is a slowly varying function. The source of the stiffness in!a stiff problem can be found only with the solution in hand. This was true for CHEM9, and again true here. It is not necessary to perform this analysis to solve a stiff problem. However, the knowledge thus gained can be very useful in understanding the nature of a specific problem. On a particular occasion, we were given a stiff problem for which the investigators expected an oscillatory solution. We determined after some analysis that the variable, which was the source of the stiffness, was the very same variable expected to change sign. This c.ras a clear indication that oscillations were impossible in that problem. If this variable were to change sign, the solution would become a positive exponential that would increase without bound. If the negative exponential in our model stiff problem is made positive, the problem is not stiff. It simply blows up. Given in TABLE 1, k2=2,85x10g for 600K; therefore, with the initial coefficient is ~=0.285. condition 12(O)= iO-lo, the expected stiff exponential This is in perfect agreement with the calculated x given in TABLE 2. For 800% k2=3.89x10g given in TABLE 1, and x=0.389 in TABLE 2. The Barnes problem k2 =4.67x10’ given in TABLE 1, and is not stiff at lOOOK, observe that x=0.47/0.7 in TABLE 2.
264
We have analyzed all stiff problems known to us, using the analysis above. In all cases, we have succeeded in obtaining agreement between the expected stiffness x and its calculated value. It is quite easy to identify the source of the stiffness in a problem using the ATOMFT system. Solving Stiff Problems Since the coefficient x is nearly constant in stiff problems, the Taylor series method is simply to factor out the exponential function by subtracting a portion What of the Taylor series terms attributable to the negative exponential. remains will be a polynomial, which is summed for analytic continuation. The numerically unstable portion of the Taylor is removed from the solution process. The manner:-
Taylor
series
method
deals
with stiff
problems
in the following
1. The length of the Taylor series used in the solution is reduced from the usual 30 terms to only 15 terms, or less. The reason for this is the numerical instability in stiff problems. 2. For every component of the stiff problem. the value of x is calculated using the approximations given in Eqs.[6] and [7]. 3. If the value of A is not constant with respect to the order n, the usual radius of convergence analysis is performed and the normal integration step is taken. For this integration step, the problem is not stiff; the rest of this algorithm is then ignored. 4. If x is found to be constant from the 15-th term down to the N-th term of the Taylor series, that series is truncated to N-l terms. We do not use the portion of the series from N to 15, where there is numerical instability. 5. The exponential function, exp(-At), is subtracted from the truuncated Taylor series starting at the very first term. Each variable may have a different x and a different N. The stiff solution is carried forward to the following steps. 6. The residual series terms are examined to make a conservative estimate for an allowable stepsize. This estimate can be many orders of magnitudes larger than the stepsize that would result from the full (untruncated) series. 7 The allowable stepsizes of all the solution variables of the stiff problem are compared, and the smallest one is used to determine the actual integration stepsize. 8. Using the optimum stepsize, analytic continuation is applied to each of the solution variables. Also, the exponential functions, exp(-At), that have been subtracted in step 5 are added back. This completes the integration step.
265
tested have been solved successfully by this algorithm. In the case of some problems, the complete so!ution is obtained with less than 25 integration steps. One of the Enright-Hull set of stiff problems, E2, was solved by ATOMFT in one step?
All the stiff problems
Stiff-Orders
Greater Than One
As defined earlier, the stiff-order of a problem is the order of the hidden lipear ODE, which generates the ‘stiff’ exponential function. Therefore, an example stiff problem with stiff-order two can be created by adding function ODE to the first Painleve transcendent equations as below.
!% dt
=
d2S iii?’
$f=
p*
-AS
+
6g2
+
t
+
a sine
s.
C$,
with g(U)=s(O)= 1. and p(O)=0 on [U, I]. With various values for x 2nd c, the problem in Eq.[ 101 can be made stiff. Analysis of the Tayio: series for problems with stiff-order two differs from Eq.[b] only with a small change in the recursion relation:
gl.papL_I.
[Ill
n
This idea can he carried forward la problems with stiff-order three and to probiems with stiff-order four. We have tested these classes of stiff problems are using examples created as above. In every instance, the ‘stiff’ coefficients found in complete agreement with theory. An ex&l;+le stiff system of ODE’s for stiff-order
with g(O)=s(O)= four is:Fl;‘(ni) n-
four is:-
I. and p(O)=0 on [O. I]. The recursion
*-l)(n-2)(n-3)(n-4) h’
= _ ). .
The case of stiff-order four is of qome ;ntcre4t because exhibits this characteristic.
relation
for stiff-order
1121
the Krogh prohlem(l*)
266
Fred Krogh”‘) very ingeniously created an artificially stiff problem by combining four first-order ODE’s in a rather complicated convolution. This problem is very interesting to analyze. It can be stated as follows:dc! =
- t- Pl
dt 2
+
p2
+
p3
=
-
(PI
-
dv\ * dt
-
(PI
+ P2
-
(PI
+ P2 + P3 -
!-ix=‘
dt
with the associated
P‘V2
+
p2 + P3 + P&W
-
P3 + P‘W
variab!es
P.e)/2
=
(v,
_
+
v2 _
-
+ g2
*
+ g3
*
+ R4
*
given by
=1 =(-yi+Yz+VJ+~*)/2 23
+ E.1 .
*
Z2=().1‘Y2+Y3+Ylb)/2
P
p3 + ?r,)/2 e 24 = (y-1 + y2 + Y3 - y/J/2
9
p1 = -1000 z1 - 1000 22 . p&a= 1000 21 - 1000 22 . p3 = 10 23 . PI. = 0.001 26 . SC:= (21 21
-
22
22)/2
s2
.
s1+s2+s3+sJ/2
gr=(-
g3 = (s1 + s2 -
s3
+
=
Zl
22
.
s3
=
z3
23
.
S&
=
, g,= (s1- s,+s,+sJ/2
S&)/2 . g, = (s1 + s2 +
s3
-
26
24
.
.
sq)/2 .
and with the initial conditions y1=
0 -
y:z
=
-2
.
y3
=
-1 , y-4= -1
The solution of this problem is sought from t =0 to t =400. The solution of this problem near t=O requires very small integration steps, because there are two pairs of simple poles very close to the ositive real axis, The first pair of simple poies are at t=5.59x10-3 + i1.66xlO- ej. The second pair of poles are at t =8.73x1O-3 -+i 1.48x10e3. The solution of the Krogh problem for times from t =0 to t =0.014 is shown in Figure 4. Only the variables yl and y3 are shown, because the other variables, y2 and y4, more-or-less shaddow these two. This problem is not stiff while its solution is under the influence of pairs of singularities. In this case, the Taylor series of the solution cannot yield a constant value for x of whatever stiff order. As shown in Figure 1. the Taylor terms describe a damped sine wave while there are singularities close by. The solution of the Krogh problem rem&ins ‘under’ the second pair of simple poles until about t =0.016. Once tne solution has passed beyond t=0.016, there gradually appears a consistent value for 1 of stiff order four. When the solution has reached t =0.025, there is a definite constant value of x=-4~10~~ for stiff order four.
267
t
1 inear vertical scale
simple pole 4 0
a_
000
0
0
0
0
0
Y,(t) -
y,(t)
0
simple poie 0 0
0 0
0
Figure 4. Solution of Krogh’s Problem. This is a very large vabe for A, and normally it would make rhe problem extremely stiff and difficult to solve. However, the solution of the Krogh problem has becomed steady state in nature. The four dependent variables show definite inter-relations; y2 = yl - JOOO, and y4 = - y-3. Beyond t =O 06. the ATOMFT solution progresses without the slightest difficulty. It is not necessary to invoke the stiff algorithm for this problem, The Krogh problem is solved easily with the usual Taylor series method for non-stiff problems. The only interest in this problem is the presence of the stiff order four A. Conclusions Our analysis of long Taylor series for the solutions of stiff problems has resulted ill an understanding of their true character. We have found a linear ODE hidden within the stiff systems of ODE’s, which yields an exponential component in the solution. Most stiff problems are of this type. Artificial stiff problems have been created by combining linear ODE’s with a non-linear Our analysis of specific stiff problems has problem with known solutions. allowed us to identify the source of the stiffness in them. So far. we have yet to encounter an ODE system which has a stiff order higher than one, and which The Taylor series algorithm for requires the stiff algorithm for its solution.
solving stiff problems is an integral automatic solution of ODE’s
part
of the ATOMt;‘l’ system
for the
Acknowledgements The author would like to express his gratitude to Fred Krogh for some very inreresting dialogues. and to Robert Corless for some critical comments on this paper. He gives particular thanks to George Corliss for Joint research on all aspects of Taylor series. References 1.
R. H. Barnes. J. L. McFarling. J. F. Kircher. and C. W. Townley. “Analytical Studies of Methyl lodide Formation under Nuclear Reactor Accident Conditions,” Rep. BMI-1816. Battelle Mem. Inst.. Columbus. Ohio, ‘.967.
2. D. Barton. 1. M. Willers. and R. V. M. Zahar. “The Automatic Solution of Systems of Ordinary Differential Equations by the Method of Taylor Series,” Computer J., v. 14, 1971, ~243. 3. Y. F. Chang, “Compiler Program Solves Equations.” Computer Decisions, June 1972, p 59: and “Automatic Soiution of Differential Equations.” Lecture Notes in Math. #430, Springer-Verlag. 1974. pp. 61-94. 4. Y. F.
Blvd.
7. Y. F. Chang and G. F. Coriiss, “Ratio-1lke and recurrence relation tests for cox’ergence of series,” J. Inst. Math. Appl., v. 25, 1980, pp. 349-359. 8. G. F. Corliss and Y. F. Chang, “Solving ordinary differential equations using Taylor series,” ACM Trans. Math. Soft., vol. 8, 1982, pp. 114- 144. “Integration 9. C. F. Curt& and J. 0. Hirschfelder, Proc. Nat. Acad. Sci., 38. 1952, pp. 235-243. IO. G. Dahlquist, “A Special Stability Methods,” BIT 3, 1963, pp, 27-43
Problem
for
of Stiff Equations,”
Linear
Multi-step
II. P. H. Davis, G. F. Corliss, and G. S. Krenz, “Bibliography on Methods and Techniques in Differentiation Arithmetic,” Technical Report No. AM-88-09, School of Mathematics, University of Bristol, 1988.
269
12. W. H. Enright, T. E. Hull, and B. Lindberg, “Comparing Methods for Stiff Systems of ODE’s, BIT 15,1975, pp. 10-49.
Numerical
13. W. H. Enright and T. E. Hulk, “Comparing Numerical Methods for the Solution of Stiff Systems of ODE’s Arising in Chemistry,” in Numerical Methods for Differential Systems, L. Lapidus and W. E. Schiesser (Eds), Academic>ess, New York, 1976, pp. 45-66. 14. A. Gibbons, “A Program for the Automatic Solution Equations using the Method of Taylor Series,” Computer pp. 108-111.
of Differential J., vol. 3, 1960,
15. H. J. Halin, “The Applicability of Taylor Series Methods in Simulation,” in Proceedings of the 1983 Summer Computer Simulation Conference, Vancouver, B. C., 1983. 16. P. Henrici, Applied and Computational New York, 1974.
Complex Analysis, Vol I, Wiley,
17. H. Kagiwada, R. Kalaba, N. Rasakhoo, and K. Spingarn, Numerical Derivatives and Nonlinear Analysis, Plenum Press, New York, 1986. t8
F T. Kroeh. “On Testing a Subroutine for the Numerical Integration Ordinary Dltferentlal equations,” J. ACM 20,1973, pp. 945-562.
19. J. A. Leavitt, “Methods and Applications vol. 20,1966, pp. 46-52. 20. R. E. Moore, 1966.
Interval
Analysis,
of Power Series,” Math. Comp.,
Prentice-Hall,
Englewood
21. R. E. Moore, Mathematical Elements of Scientific Rinehart and Winston, New York, 1975. 22. R. E. Moore, 1985.
Computational
Functional
23. L. B. Rail, Automatic Differentiation: Springer, New York, 1981.
of
Analysis, Techniques
Cliffs, N.J.,
Computing,
Wiley,
Holt,
New York,
and Applications,
24. L. F. Shampine, “Evaluation of a Test Set for Stiff ODE Solvers,” ACM Trans. Math. Soft., vol. 7, 1981, pp. 409-420. 25. R. E. Wengert, “A Simple Automafic Comm. ACM, vol. 9, 1964,443-464.
Derivative
Evaluation
Program”,
Solving STIPF systems
by Taylor
Series
Y. F. C-hang 976 U
Foothill BLvd.. Claremont. CA. 91711
Abstract The analysts ot long Taylor series for the solutions of stiff prohlems gives new insight into their true character and leads to new solution algorithms. Some ODE’s are difficult to solve becau se of singularities close to the real axis. This is not true for stift problems. There is a linear ODE hidden within the stiff system of ODE’s. When the hioden linear ODE is of order one, there is an exponential component in the solution. Most stiff problems arc of this type. This exponential forces the integration stepsize to be very small relative to the domarn of the stiff problems. This exponential component is readily made apparent by analysis of the Taylor series for the solutions. We analyze in detail a stiff problem from chemistry, which is very difficult to solve. The portions of the equations giving rise to the stiffness in this problem are identified. In this problem. there are more than one hidden linear ODE. Artificial stiff problems with up to stiff-order four have been created by combining linear ODE’s with a non-linear problem with known solutions. The Taylor series algorithm for solving stiff problems is an integral part of the ATOMFT system for the automatic solution of ODE’s. Introduction:-
Taylor
Series Methods
The use of Taylor series to solve ODE’s began with Sir Issac Newton, who of solved ODE’s with 4-term Taylor series. In modern times, the development Taylor series methods for ODE’s was independently pursued by Gibbons(14), Moore(21), Leavitt(“), Barton, et al(l), Chang, et al(3*8), Rail(“), and Halin(15). They were all able to demonstrate that Taylor series methods do make a practical algorithm for the solution of ODE’s. This author’s Taylor was shown to be faster than the usual series package called ATOMCC(*) methods when solving the problems in the Toronto non-stiff ODE test set at high precisions. There is at leact one book (**I that advocates the Taylor series methods as general purpose algorithms for ODE’s. There are at least four misconceptions concerning the efficiencv. stability, and general. viability of the Taylor series methods as applied to ODZ’s. They ale:- (1) difficulty to evaluate the high order derivatives, (2) sequential differentiation by the use of recursion relations is grossly unstable, (3) Henrrci’s ‘proof’(16), and (4) a ‘proof’ that the high order derivatives may be unbounded. The facts speak!: against these misconceptions. Analyses of the high order Taylor series terms in ODE problems yield the precise location and the precise @ELsevier Science Pubhshing Co., Inc., 1989 655 Avenue of the Ameruxs, New York, NY
10010
0096-3003/99/$03.50
252
order of the singularities in the solution functions. Such precision would be impossible without numerical stability. We shall answer these misconceptions in the following paragraphs. ( 1)Many numerical analysis texts echoes the incorrect conventional wisdom that Taylor series methods are impractical due to difficulties in computing the Longs Tay&r2;e;iy are sequentially evaluated using high order terms. with accuracy and efficiency. It is not differentiationarithmetics - * a svmbolic process. and it does not use finite differences. There is a recent hibiiography(ll) containing 90 references to differentiation arithmetics and its applications. The recursion relations used in differentiation arithmetics are derived from the OPE’s themselves by applications of the chain rule. There is no difficulty in evaluating high order Taylor series terms. l
(2)There are those who hold to the view that the sequential evaluation of Taj-lor series terms is grossly unstable and leads to unusable series terms. Although it is true that all sequential processes are inherently unstable, this is not of any major concern in the evaluation of Taylor series terms, because the magnitude of the instability is usually very small. Except for one class of problems. w-e have found ail recursion relations to yield very good results. A Taylor series of 10.000 terms have been calculated with high precision for the eIliptic-integral function. The one class of problems, where the instability does cause some concern, are called ‘last-pole’ problemsC5). This error condition occurs when the solution in these problems approaches an aysmptotic entire function, which does not have any singularities in the finite complex plane. The difficulty arises as the solution passes beyond the ‘last-pole’ in the non-linear domain of the problem into the linear domain. There are subtractive errors that arise from the ODE’s becoming almost linear with a vanishingly small (and so the subtractive error) non-linear component. (3)It was unfortur.ate for Henrici(16) to have offered a proof that analytic continuation obtained by rearranging a truncated Taylor series is unstable. Although his proof is correct, it is totally irrelevant. None of the Taylor series methods developed for ODE’s performs analytic continuation by rearranging the terms of a truncated series. They ail generate a new Taylor series for every function at each new integration step. They all have successfully solved ODE’s through hundreds of integration steps extending far beyond the domain of Henrici’s concern. (4)There is another ‘proof’ that is doing the rounds of numerical analysis texts. It deals with the fact that high order derivatives increase without bound when there are singularities close by. This authorC4) has shown this proof to be correct but also totally irrelevant. The derivatives do increase without bound; however. the Tqlo: series terms are bounded. As a matter of fact, the Taylor series terms typically describe a damped sine wave. Figure 1 is a graph of the first 100 terms of F(n) versus n at x=0 for the solution off” + f’f + f = 0, withf(O)=-2.81 andf’(O)=-0.99. There are simple poles at 2=2.2 f i0.69; the stepsize used (h=2.3) is slightly less than the radius of convergence. The angle of the poles measured from the real axis at x=0 (17.5”). is exactI) equal to the phase angle difference between the Taylor series terms.
Figure 1. Typical Variation
of Taylor Series Terms from ODE’s
The Taylor series methods implemented by each of the Barton, Chang. and Halin teams were compilers that accepted mathematical statements of the differential equations to be solved. The compilers produced an object program that is executed to solve a specific problem. Thus. problems are solved ‘tiatomatically’: no formal programming is necessary. This represents a significant advantage over the usual methods. The original compiler(3) written by this author was computer specific. it Further was superceeded by a portable program called ATOMCC(4*s’. extensions to include the solution of stiff problems, arbitrary functions definable by the user, and the soiution of DAE’s have led to the present ATOMFTc6) program. The Taylor series method 4as evolved into a general purpose numerical technique capable of solving all manner of problems. The ATOMFT packagec6) applies the detailed analyses(3*7) of the Taylor series terms to advance the understanding and treatment of the truncation Tne Taylor series analyses yield the location errors in numerical integration. and order of the singularities in the solution functions. Thus, ATOMFT is able to control local errors with sufficient accuracy to achieve proportionally between the local trunation and the global errors. This also represents a significant advantage over the usual numerical integration techniques. Introduction:-
Taylor Series Solution of Stiff Problems
Stiff s stems of ODE’s were first discussed by Curtiss et alc9’ in 1952. as the source of the difficulty Dahlquist(“) had identified numerical instabilit conventional Enright, et al (12.13 examined in solving these problems. numerical methods for stiff problems. The Taylor series method to be discussed is different from conventional methods. First, we use 30-term Taylor series; thus, it is a 29th order method. Second, we analyze the series to obtain the radius of convergence and other characteristic constants such as exponential coefficients. Third, the integration stepsizes used are much larger than in other methods and they are infinitely variable.
254
In 1967, while developing the analysis of Tayior series, we were given a problem by Barnes, et al CL). This problem was stiff. It did not obey any of the Recent studies has revealed simple rules or formulae being then developed. recursion relationships identifiable with stiff problems. The numerical instability inherenr in the recursively generated Taylor series terms does I.ot show up in any casual inspection of the Taylor terms. In contrast to the last-pole problems mentioned above, most if not all non-linear problems have Taylor series terms that are reliably reproduced. Rather than the propagation of large subtractive errors as in last-pole problems, the numerical instability in stiff problems is quite difficult to detect. We need to ha-fe a clear understanding of this instability in order to accurately and reliably evaluate an exponential coefficient A. The singificance of this exponential coefficient is explained in our model stiff problem in the next section. In our solution algorithm, once the coefficient A is found to be constant, we truncate the Taylor series to remove ALL of the terms that are numerically unstable. Therefore, we have the best of both worlds. The numerical instability is analyzed to find a constant coefficient A. Then, that value of x is used to take giant integration steps, while all those Taylor terms used to evaluate that x are removed. A Model Stiff Problem
When long Taylor series for every function in a problem is available for analysis, simple relationships between adjacent Taylor terms become readily apparent. Assume, for the moment, that a stiff Frcblem contains a first order linear ODE within a system of non-linear ODE%. This assumption will be justified by the analysis of data from actual stiff problems. Therefore, we have
whose solution is the exponential function, s(t)=exp(-At). the total solution for the stiff problem is of the form
f(t) wheref(t)
=
S(t)
+
g(t)
We also assume that
,
is the total solution, and g(t) arises from the non-linear
[21 stiff ODE%.
The adjacent Taylor terms of the S-series are related by:-
*
(
n-
l)e+=-x,
where n is the order of the terms, and h is the stepsize. When the magnitude of the exponential coefficient, x, is large, the magnitude of S(n) increases very rapidly with n. The magnitude of S(n) reaches a maximum at some very large value of n, and then decreases for all higher n’s, The peak of this ‘hump’ for the
usual stiff problem occurs beyond beyond our sphere of interest.
hundreds
or even thousands
of terms,
far
The magnitude of the Taylor term, G(n), is determined by the positions of the singularities in the solution in the complex plane of the independent variable t. The radius of convergence is defined to be the distance from the point of expansion to the nearest singularity. If the stepsize is equal in magnitude to or less than) the radius of convergence, the Taylor term, G(n), is bounded(5* s 1. Typically, the Taylor term, G(n), describes a damped sine wave (see Figure 1). This is true for solutions of -all non-linear ODE systems. Since the Taylor series for g(t) is bounded. and the function s(t) has terms that increase very rapidly with n, the high order F(n) are dominated by S(n). Therefore, if x is large, we cau expect the same relationship between adjacent Taylor terms of f(t) as for s(t) when n is large, say around IO.
P- 1)y+--a (
n-
When A is large, the ratio F(n)/F(n-1) is large except when the stepsize is very small. This is the source of the stiffness in stiff problems, based on the assumption stated in Eqs.[ I] and [2]. When solving a non-stiff non-linear ODE system, the stepsize is restricted to be less than the radius of convergence. If there are singularities close to the real axis, the solution must proceed with small stepsizes. Any attempt to take large steps will result in solutions that have very large errors. In contrast, the solution of a stiff system of ODE’s has a stepsize forced to be small by an exponential component. In this case, taking giant steps should not result in appreciable errors if we account properly tnr the singularities in the non-linear g(t) function. This is our model of stiff prol!tms. We will attempt to fit this model to some actual stiff problems. The stiff-order of a stiff problem is defined to correspond to the order of the hidden linear ODE. Almost all of the non-linear problems in the EnrightHull set(13* 14) with and without the modifications recommended by Shampine(17) have stiff-order one. We will examine stiff problems with stifforder higher than one after completing the analysis to find the exponential coefficient x.
Detection
of Exponentials,
Stiff-Order
One
In this section, the analysis of Taylor series problems in the search for a reliable evaluation of the The numerical instability discussed by Dahlquist(“) We have found it necessary to analyze this instability evaluate the exponential coefficient x.
is applied to many stiff exponential coefficient X. is very much in evidence. in full in order to reliably
Barnes(l) problem was concerned with the production, inside a run-away nuclear reactor, of radioactive methyl-iodide (CH3I) from methane (CH4), free
iodine (12). hydrogen iodide (HI). the iodine ion (I). and the methyl ion (CH3). The iodide reaction equations for the nuclear accident problem of Barnes are:n I2 _-.
dt
- - kj 12 # + k6 1' H .
d CH4
- dt
d CH3I -__--dt
_ k3 CH3 HI - k, CH4 I ,
141
= _ k: CH31 I + kZ CH3 I2 .
M = I + CH4 + CH3 + HI + 12 + CH31 .
with the initial concentrations
given by
2 12(O) = 2 12 + HI + I + CH31 . .G CH4(0) = 4 CH4 + 3 CH3 + HI + 3 CH3I . and CH4(0) = CH4 + CH3 + CH31 .
The initial concentrations
are:- CH4(0)=4x10‘8.
The reaction rate constants. k’s, depend nuclear reactor. They are given in Table 1. TABLE Temperature
%% 600K
kl
8.31E+6 6_30E+5 8_7OE+3
1. The Reaction kz
4.67E+9 3.89E+9 2.85E+9
12(O)= 10‘lo, and CH31(0)=0. on the temperature
Rate Constants. ks
ks l_OOE+9 7.43E+8 4.57&+8
inside the
2.46E+4 3.55E+2 0.317
kg
5.51E+4 1.43E+9 ;.9887v2 1.92E+9 2.84E+9 .
The numerical solution of this problem at the temperature of IOOOKis quite easy. The problem is hardly stiff at all. and it is easily solved by most numerical methods. At the temperature of 80OK, this probiem is stiff but not difficult to solve. At the temperature of 6OOK. this problem is extremely stiff. The chemical reaction times to reach steady state inside the nuclear reactor are:- 10,000 seconds for LOOOK,300.000 seconds for 80OK, and 60,000,OOO seconds for 600K. These are chemical reaction times and NOT computer calculation times. Under the assumed model stiff problem of Eqs.[ I] and [2]. the ratio given in Eq.[3] should yield a coefficient A that is constant. We analyzed the Taylor series for 12 in the three temperatures for the Barnes problem. Listed in TABLE 2 are values of x for all Taylor orders from n=7 to n=30 at a fixed time. These ratios are calculated with double-precision; therefore, they should be accurate to at least 12-digits. For our proposed model of stiff problems to be valid, the
257
numbers listed in TABLE 2 should be constant the Taylor order, n. TABLE 2. The Exponential n
lOOOK :
BOOK
.389177 “4:29:3: .389214 .389250
:%26d: .389286 .485546 .487678 .489805 491926
: %Ei
:"499% .498290 .500440 .5026?9
.I%9500 :388%%
: 5500% .510336
-389616 .389623 .389601
.513975 .519020
.526669 .539000 .559320 .592065 :sM;
: %k% .389465
.389595
: %53% .388892 .387993 : %4% .374558 .357958
in each column independent
Coefficients, 600K
of
x. Krogh
1000.03 999.99 999.99 999.99 Ext 999:92 999.84 999.69 999.38 E% ;;;:cI'j 979:a7 958.93 914.36 812.71 539.24 -5E ii 2654.85 2244.88 2106.70
In fact, the x’s calculated are not ctinstant to the degree expected, casting doubts on the validity of the model. Further study is needed to establish the truth of this model. Also, if this model is correct, the stiffness is not determined by the magnitude of the exponential coefficient, x. Most of the coefficients in TABLE 2 are smaller than 0.5. The 6OOK problem is the most stiff among the three, and yet it has the smallest coefficient. Is the stiffness of a problem indicated by the constancy of the coefficient. 3 This is shown to be untrue by the results [in the last column of TABLE 21 for the very stiff Krogh problem(18) (also E4(12)). The calculated coefficients for Krogh’s problem are certainly not constant. Dahlquist(l’) had identified numerical instability as the source of the it behooves us to suppose that difficulty in solving these problems. Therefore, the non-constaucy of x shown in TABLE 2 is due to an accumulation of errors. For the moment, let this error propagation be linearly dependent on the order of the terms. We will fit the data by a straight line x,
=
x,
+
an,
161
where X, is the ratio listed in TABLE 2, x0 is this ratio extrapolation to the zeroth term (the vertical intercept), n is the order of the Taylor term, and o is the slope. Results of applying this fit to the Barnes problem and many other stiff problems indicate that the numerical instability is more complex than a simple linear propagation of errors.
258
We next fit the non-constant x in TABLE 2 to a propagation of errors that is plus an also dependent on the power of n. This fit is a linear correction exponential correction; &I
=
.b
+
an
+ ar”
i7]
.
where 6 is a constant, and $’ is the exponential component With this approximaticn, all the stiff problems examined constant with $-digit precision.
of the correction. yield A’\ that are
The results for the Barnes problem at 80OK are listed in TABLE 3. The first three terms cannot be fitted because the linear exponential function. s(t) has yet to be dominant. The next six terms have a linear behavior. All the remaining terms required the full approximation. These results are typical of all the stiff problems examined. TABLE 3. Calculated n 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
raw X
error
-3.8917D-03 4.3D-02 -3.8921D-01 9_6D-05 -3.89250-01 9.2D-05 -3.89280-01 1.8D-07 -3.8932&01 -3.89351)-01 -3.89398-01 -3.8942D-01 -3_8946D-01 -3.895OD-01 -3.89531)-03j-g:;; -3.8956D-01 7:aD-oa -3.8959D-01 8.4D-08 -3.a961D-01 -3.89621)-01P*s% -3.8960D-01 2:5D107 -3.8952D-01 8.2D-07 -3.8932D-01 3.38-06
Exponential fit X
Coefficients.
Q
-3.88998D-01 -3.60D-05 -3_88998D-01 -3.59D-05 -3_8899aD-01 -3.59D-05 -3.889991)-01-3.598-05 -3.89000D-01 -3.57D-05 -3_89003D-01 -3.55D-05 -3.88998D-01 -3,60D-05 -3_88998D-01 -3.6OD-05 -3.88998D-01 -3_6OD-05 -3.88998D-01 -3.60"-05 -3.88998D-01 -3.608-05 -3.88998D-01 -3_60D-05 -3.88998D-01 ->.60D-05 -3.8899aD-01 -3.59D-05 -3.89000D-01 -3.59D-05
80OK.
B
g:;:
-7
z%z:
17D-10 z:oD+oo 15D-10 2,oD+oo z.oD+Oo gr,l:: 2.oD+oo 13D-10 12D-10 EE% llD-10 2:oD+oo
The important conclusion of this study is not our ability to the fit the data with a complex error analysis; it is the clear existence of a constant A. In the fourth coIumn, the value of A &es by 5 in the sixth digit. Therefore, we conclude that the model for stiff problems, Eqs.[l] and [2J, is valid. A problem is stiff if upon application of Eqs.[6] and ]7] a constant A is found. Incidentally, the accuracy of x does not materially effect the accuracy of the ODE solution. As is explained later in the section on the stiff algorithm, the exponentia1 function, exp(-at), is first subtracted from the solution and then added back at a later stage. It bears repeating that, in our solution algorithm, once the coefficient x is found to be constant, we truncate the Taylor series to remove ALL of the terms that are contaminated by numerical instability
259
Synthesis
of a Stiff Problem
In this section, we further demonstrate the validity of our model for stiff problems by synthesizing a stiff problem. The simplest model of a stiff problem consists of a system of non-linear ODE’s for the function g(t) plus a linear ODE, Eq.[ 1). Hence, it should be possible to create a stiff problem artificiallv. Until we can create a stiff problem out of known components, we cannot claim to understand stiff problems. This task is undertaken equations can be as follows:-
c&I
dt
-
p
s!.E ’
dt
=
in small steps.
6 g”
+
t
+
s ,
First, a possible
stiff system of
and
with g(0) = s(0) = 1, and p(O) =0 on [0, 11. The solution to the first two equations, Eq.[8] without the s at the end, is the first Painleve transcendent. To this well known problem, we add the function, s, and the associated linear ODE. The first Painleve transcendent, g, is a monotonicaly increasing function with a simple pole on the real axis at ? = 1.2. The addition of the linear function, s, to the Painleve transcendent does not As in any non-stiff solution, s(t), rapidly make this system a stiff problem. decays to some negligible magnitude for a large x. Therefore, beyond a very small step in time, the contribution of s(t) to the solution g(t) becomes negligible, and the solution proceeds ignoring the exponential. Therefore, this simple combination of ODE’s does not produce a stiff problem. In a stiff problem, as defined by our model, insignificance. It must somehow persist in disturbing the in spite of its rapid decay. Therefore, there must be a the non-linear ODE’s for g(t) and the linear ODE magnitude for s(t) can remain large. To accomplish replace the simple linear ODE by ds dt
= - A s
+
c g(t)
,
s(t) cannot decay to solution for large times cross-coupling between for s(t) such that the this cross-coupling, we
191
where c is the cross-coupling constact. This system of ODE’s is indeed stiff for non-vanishing values of c. We have succeeded in creating a stiff problem from known components. The stiffness of this problem is controlled by the magnitude of the coefficient, x, and the amount of the cross-coupling. After studying the results from this synthesis model, we conclude that the magnitude of x does not by itself determine the stiffness. Some very stiff problems have x’s that are quite small. The magnitude of the cross-couplmg constant isvery important in the stiffness of a problem.
260
Analysis of Stiff Problems In this sectioo, we take the final step towards a clear understanding of stiff problems. Since the model for stiff problems contains a linear ODE which has a cross-coupling component from the non-linear ODE’s, it should be possible to ideniify and isolate the source of stiffness in problems. It is important in this identification to remove all non-linear portions of each ODE to see whether it satisfies the form of Eq.191. The coefficient A is allowed to be a slowly-varying function of time. The final test, to identify the source of the stiffness in the problem, is to compare the value of the slowly varying function against the value of a calculated from the Taylor series. Of al! the stiff problems known to us, the most difficult ones are the Barnes methyl-iodide problem and CHEMP in Enright-Hu11(13). We now identify the sources of the stiffrzss in CHEMO. This problem has two distinct stiff phases, separated by phases of uormal behavior. The ODE’s for CHEMP are:-
ax
dt
= 77.27(y - xy + x - 8.3758-6 2)
xy + z)/77.27 ,
Q=(_y_
df ZL
,
and
= 0.161(x - 2) .
with initialvalues
x(O)-4,y(O)=
1.1, r(0)=4.
on [0,300].
To identify the source(s) of a problem’s stiffness, we look for occurrences of linear ODE’s. It is important to allow the exponential coefficient, A, to be a function of t. Only after having solved the problem, and having evaluated that function. can we determine if that function is nearly constant. The procedure is to force the ODE to fit the following form, collecting all the relevant portions into the X(t) and all the irrelevant portions into the [whatever]. ds zG=
- A(t) s
+
[whatever]
With three ODE’s in CHEM9, there are three potential model. They are:(1)
dx ;ii:=
(2)
d-1!= dt
-
(3)
dr z=
- (0.161) z
After
candidate,
- 77.27 (y - 1) X
$?$.
y
+ f
+
occurrences
of the above
[whatever] ,
[whatever] .
and
[whatever]
having solved CHEMY, we can eliminate ease(3) as a potential because at no time in the solution is there a x=0.161. That leaves
261
both ca\e( I) and cu\e(2) a\ possihle contenders for the source of the stiffness. Case( I ) has x(t) = 77.27 (v- I ). Therefore. .v(t) must he slowly varying for case( 1) to he the source. Case(Z) ha> x(t)=(x+ 1)/77.27. Therefore, x(t) must he slowly varytng for CHEMO to he stiff. So. the search is for phases of the solution when either u(t) or l:(t) I\ varying slowly. The complete
of CHEMV i\ shown in Figure 2. There is a sharp The vartation of y(t) is somewhat more gradual. The two functions are shown together in Figure 2: however, they differ in magnitude hy a factor of 70. The spike in x(t) is 70 times larger than the peak of r(t). An enlarged view of the solution near the origin is shown in Figure 3. spike
in .t(t)
solutton
near the origin.
linear -- vertical scale *. 9 .
.
oocl
.
.
-
x(t)
-
Y(t)
-. -. -. -... --.__ . ...**. 0
100
0
Figure 2. Solution of CHEMY. for x and y.
There are four distinct phases in the solution of this problem. The four phases are in the times [0, 1.256). (1.256.3 g), (3.8.3.98), and (3.98.300). 1. In the first phase IO. 1.256), our Taylor series solution shows a conjugate pair of singularittes at z=-0.35 + 10.30 and another pair at z= 1. I5 + iO.03. The second singularity pair causes x(t) to increase exponentially from x(0)=4.0 to x( 1.256)= 1.1784~10~. Here, CHEM9 is not a stiff problem. Whenever there is a rapidly varying function in a problem, that problem cannot he stiff. This is because a rapidly varying function is a clear indication that the solution is under the influence of some singularities. Then, the Taylor series of the solution cannot yield a constant value for X. As shown in Figure 1. the Taylor terms describe a damped sine wave when there are singularities close by.
262
In the second phase (1.256.3.8) x(t) is slowly varying. Therefore, ease(2) is the source of the stiffness. The calculated A agrees with the predicted value of (x+ 1)/77.27. The magnitude of x(t) decreases linearly from x(1.256)=1.1784x105with~=1.500tox(3.8)=4.968x104withx=692. 2. In the third phase (3.8, 3.98), x(t) decreases very rapidly from x(3.8)=4.968x10* to x(3.98)= 1.04. CHEM9 is again not a stiff problem. 3. In the fourth phase (3.98.3001, y(t) is slowly varying. Therefore. case( 1) is the source of the stiffness. The calculated x agrees with the expected value of 77.27 (y-l). The magnitude of y(t) increases slowly from y(3.98)=29 with x=2.2x103 to a maximum of y(17.1)=1.757xld with x = 1.36x105, and then decreases slowly to ~(300) = 1.22 with x = 17.
t
linear vertical scale stiff Case(P) stiff Case(l)
000
-
x(t)
. . . - Y(t)
5
6
7
a
9
10
Figure 3. Solution of CHEM9, an Enlarged View. The regions of stiff solutions are marked in Figure 3. It is clear in this graph that the spike in x(t) has a definite structure. The first stiff region does not begin until after x(t) has reached its peak and is decreasing on a sloped plateau.
263
The first stiff region ends when x(t) drops too rapidly !o sustain a ‘constant” x for Case(Z). The second stiff region begins once y(t) has becomed large enough to make the x for Case( 1) a dominant factor. The stiff algorithm is ahsolutly necessary in the solution of CHEMY. In the first stiff region. the integration stepsize!, are about ten times larger than a normal Taylor series stepsize. In the second stiff region. the integration stepsizes are at least a thousand times larger than the normal stepsizes. We now attempt to identify the stiffness in Barnes problem. Refering to the ODE’s in Eq.[4], if 12 is to be the linear function. then (M) must be slowly varying, but it is not. If CH4 is to be the linear function. then (HI) must be slowly varying. but it is not. If CH3I is to be the linear function. then (I) must be slowly varying, but it is not. Therefore, we must look elsewhere. After considerable following substitutions.
search, we found CH31 to be the linear function From Eq.[SJ, we identify CH3(t) to be
by the
CH3(t) = CH4(0) - CH4(t) - CH3I(t) .
Substituting this result into the term (k2 CH3 t2) in Eq.[7] and isolating the relevant and irrelevant parts, we find the linear first-order ODE that is the stiffness of this problem. d CH31 -.X dt
- (k,I2)CH31 +
[whatever],
where A(t)=k212. During the extensive SPYC~. we needed the availability of the solution, where it was determined that I2(t) is a slowly varying function. The source of the stiffness in!a stiff problem can be found only with the solution in hand. This was true for CHEM9, and again true here. It is not necessary to perform this analysis to solve a stiff problem. However, the knowledge thus gained can be very useful in understanding the nature of a specific problem. On a particular occasion, we were given a stiff problem for which the investigators expected an oscillatory solution. We determined after some analysis that the variable, which was the source of the stiffness, was the very same variable expected to change sign. This c.ras a clear indication that oscillations were impossible in that problem. If this variable were to change sign, the solution would become a positive exponential that would increase without bound. If the negative exponential in our model stiff problem is made positive, the problem is not stiff. It simply blows up. Given in TABLE 1, k2=2,85x10g for 600K; therefore, with the initial coefficient is ~=0.285. condition 12(O)= iO-lo, the expected stiff exponential This is in perfect agreement with the calculated x given in TABLE 2. For 800% k2=3.89x10g given in TABLE 1, and x=0.389 in TABLE 2. The Barnes problem k2 =4.67x10’ given in TABLE 1, and is not stiff at lOOOK, observe that x=0.47/0.7 in TABLE 2.
264
We have analyzed all stiff problems known to us, using the analysis above. In all cases, we have succeeded in obtaining agreement between the expected stiffness x and its calculated value. It is quite easy to identify the source of the stiffness in a problem using the ATOMFT system. Solving Stiff Problems Since the coefficient x is nearly constant in stiff problems, the Taylor series method is simply to factor out the exponential function by subtracting a portion What of the Taylor series terms attributable to the negative exponential. remains will be a polynomial, which is summed for analytic continuation. The numerically unstable portion of the Taylor is removed from the solution process. The manner:-
Taylor
series
method
deals
with stiff
problems
in the following
1. The length of the Taylor series used in the solution is reduced from the usual 30 terms to only 15 terms, or less. The reason for this is the numerical instability in stiff problems. 2. For every component of the stiff problem. the value of x is calculated using the approximations given in Eqs.[6] and [7]. 3. If the value of A is not constant with respect to the order n, the usual radius of convergence analysis is performed and the normal integration step is taken. For this integration step, the problem is not stiff; the rest of this algorithm is then ignored. 4. If x is found to be constant from the 15-th term down to the N-th term of the Taylor series, that series is truncated to N-l terms. We do not use the portion of the series from N to 15, where there is numerical instability. 5. The exponential function, exp(-At), is subtracted from the truuncated Taylor series starting at the very first term. Each variable may have a different x and a different N. The stiff solution is carried forward to the following steps. 6. The residual series terms are examined to make a conservative estimate for an allowable stepsize. This estimate can be many orders of magnitudes larger than the stepsize that would result from the full (untruncated) series. 7 The allowable stepsizes of all the solution variables of the stiff problem are compared, and the smallest one is used to determine the actual integration stepsize. 8. Using the optimum stepsize, analytic continuation is applied to each of the solution variables. Also, the exponential functions, exp(-At), that have been subtracted in step 5 are added back. This completes the integration step.
265
tested have been solved successfully by this algorithm. In the case of some problems, the complete so!ution is obtained with less than 25 integration steps. One of the Enright-Hull set of stiff problems, E2, was solved by ATOMFT in one step?
All the stiff problems
Stiff-Orders
Greater Than One
As defined earlier, the stiff-order of a problem is the order of the hidden lipear ODE, which generates the ‘stiff’ exponential function. Therefore, an example stiff problem with stiff-order two can be created by adding function ODE to the first Painleve transcendent equations as below.
!% dt
=
d2S iii?’
$f=
p*
-AS
+
6g2
+
t
+
a sine
s.
C$,
with g(U)=s(O)= 1. and p(O)=0 on [U, I]. With various values for x 2nd c, the problem in Eq.[ 101 can be made stiff. Analysis of the Tayio: series for problems with stiff-order two differs from Eq.[b] only with a small change in the recursion relation:
gl.papL_I.
[Ill
n
This idea can he carried forward la problems with stiff-order three and to probiems with stiff-order four. We have tested these classes of stiff problems are using examples created as above. In every instance, the ‘stiff’ coefficients found in complete agreement with theory. An ex&l;+le stiff system of ODE’s for stiff-order
with g(O)=s(O)= four is:Fl;‘(ni) n-
four is:-
I. and p(O)=0 on [O. I]. The recursion
*-l)(n-2)(n-3)(n-4) h’
= _ ). .
The case of stiff-order four is of qome ;ntcre4t because exhibits this characteristic.
relation
for stiff-order
1121
the Krogh prohlem(l*)
266
Fred Krogh”‘) very ingeniously created an artificially stiff problem by combining four first-order ODE’s in a rather complicated convolution. This problem is very interesting to analyze. It can be stated as follows:dc! =
- t- Pl
dt 2
+
p2
+
p3
=
-
(PI
-
dv\ * dt
-
(PI
+ P2
-
(PI
+ P2 + P3 -
!-ix=‘
dt
with the associated
P‘V2
+
p2 + P3 + P&W
-
P3 + P‘W
variab!es
P.e)/2
=
(v,
_
+
v2 _
-
+ g2
*
+ g3
*
+ R4
*
given by
=1 =(-yi+Yz+VJ+~*)/2 23
+ E.1 .
*
Z2=().1‘Y2+Y3+Ylb)/2
P
p3 + ?r,)/2 e 24 = (y-1 + y2 + Y3 - y/J/2
9
p1 = -1000 z1 - 1000 22 . p&a= 1000 21 - 1000 22 . p3 = 10 23 . PI. = 0.001 26 . SC:= (21 21
-
22
22)/2
s2
.
s1+s2+s3+sJ/2
gr=(-
g3 = (s1 + s2 -
s3
+
=
Zl
22
.
s3
=
z3
23
.
S&
=
, g,= (s1- s,+s,+sJ/2
S&)/2 . g, = (s1 + s2 +
s3
-
26
24
.
.
sq)/2 .
and with the initial conditions y1=
0 -
y:z
=
-2
.
y3
=
-1 , y-4= -1
The solution of this problem is sought from t =0 to t =400. The solution of this problem near t=O requires very small integration steps, because there are two pairs of simple poles very close to the ositive real axis, The first pair of simple poies are at t=5.59x10-3 + i1.66xlO- ej. The second pair of poles are at t =8.73x1O-3 -+i 1.48x10e3. The solution of the Krogh problem for times from t =0 to t =0.014 is shown in Figure 4. Only the variables yl and y3 are shown, because the other variables, y2 and y4, more-or-less shaddow these two. This problem is not stiff while its solution is under the influence of pairs of singularities. In this case, the Taylor series of the solution cannot yield a constant value for x of whatever stiff order. As shown in Figure 1. the Taylor terms describe a damped sine wave while there are singularities close by. The solution of the Krogh problem rem&ins ‘under’ the second pair of simple poles until about t =0.016. Once tne solution has passed beyond t=0.016, there gradually appears a consistent value for 1 of stiff order four. When the solution has reached t =0.025, there is a definite constant value of x=-4~10~~ for stiff order four.
267
t
1 inear vertical scale
simple pole 4 0
a_
000
0
0
0
0
0
Y,(t) -
y,(t)
0
simple poie 0 0
0 0
0
Figure 4. Solution of Krogh’s Problem. This is a very large vabe for A, and normally it would make rhe problem extremely stiff and difficult to solve. However, the solution of the Krogh problem has becomed steady state in nature. The four dependent variables show definite inter-relations; y2 = yl - JOOO, and y4 = - y-3. Beyond t =O 06. the ATOMFT solution progresses without the slightest difficulty. It is not necessary to invoke the stiff algorithm for this problem, The Krogh problem is solved easily with the usual Taylor series method for non-stiff problems. The only interest in this problem is the presence of the stiff order four A. Conclusions Our analysis of long Taylor series for the solutions of stiff problems has resulted ill an understanding of their true character. We have found a linear ODE hidden within the stiff systems of ODE’s, which yields an exponential component in the solution. Most stiff problems are of this type. Artificial stiff problems have been created by combining linear ODE’s with a non-linear Our analysis of specific stiff problems has problem with known solutions. allowed us to identify the source of the stiffness in them. So far. we have yet to encounter an ODE system which has a stiff order higher than one, and which The Taylor series algorithm for requires the stiff algorithm for its solution.
solving stiff problems is an integral automatic solution of ODE’s
part
of the ATOMt;‘l’ system
for the
Acknowledgements The author would like to express his gratitude to Fred Krogh for some very inreresting dialogues. and to Robert Corless for some critical comments on this paper. He gives particular thanks to George Corliss for Joint research on all aspects of Taylor series. References 1.
R. H. Barnes. J. L. McFarling. J. F. Kircher. and C. W. Townley. “Analytical Studies of Methyl lodide Formation under Nuclear Reactor Accident Conditions,” Rep. BMI-1816. Battelle Mem. Inst.. Columbus. Ohio, ‘.967.
2. D. Barton. 1. M. Willers. and R. V. M. Zahar. “The Automatic Solution of Systems of Ordinary Differential Equations by the Method of Taylor Series,” Computer J., v. 14, 1971, ~243. 3. Y. F. Chang, “Compiler Program Solves Equations.” Computer Decisions, June 1972, p 59: and “Automatic Soiution of Differential Equations.” Lecture Notes in Math. #430, Springer-Verlag. 1974. pp. 61-94. 4. Y. F.
Blvd.
7. Y. F. Chang and G. F. Coriiss, “Ratio-1lke and recurrence relation tests for cox’ergence of series,” J. Inst. Math. Appl., v. 25, 1980, pp. 349-359. 8. G. F. Corliss and Y. F. Chang, “Solving ordinary differential equations using Taylor series,” ACM Trans. Math. Soft., vol. 8, 1982, pp. 114- 144. “Integration 9. C. F. Curt& and J. 0. Hirschfelder, Proc. Nat. Acad. Sci., 38. 1952, pp. 235-243. IO. G. Dahlquist, “A Special Stability Methods,” BIT 3, 1963, pp, 27-43
Problem
for
of Stiff Equations,”
Linear
Multi-step
II. P. H. Davis, G. F. Corliss, and G. S. Krenz, “Bibliography on Methods and Techniques in Differentiation Arithmetic,” Technical Report No. AM-88-09, School of Mathematics, University of Bristol, 1988.
269
12. W. H. Enright, T. E. Hull, and B. Lindberg, “Comparing Methods for Stiff Systems of ODE’s, BIT 15,1975, pp. 10-49.
Numerical
13. W. H. Enright and T. E. Hulk, “Comparing Numerical Methods for the Solution of Stiff Systems of ODE’s Arising in Chemistry,” in Numerical Methods for Differential Systems, L. Lapidus and W. E. Schiesser (Eds), Academic>ess, New York, 1976, pp. 45-66. 14. A. Gibbons, “A Program for the Automatic Solution Equations using the Method of Taylor Series,” Computer pp. 108-111.
of Differential J., vol. 3, 1960,
15. H. J. Halin, “The Applicability of Taylor Series Methods in Simulation,” in Proceedings of the 1983 Summer Computer Simulation Conference, Vancouver, B. C., 1983. 16. P. Henrici, Applied and Computational New York, 1974.
Complex Analysis, Vol I, Wiley,
17. H. Kagiwada, R. Kalaba, N. Rasakhoo, and K. Spingarn, Numerical Derivatives and Nonlinear Analysis, Plenum Press, New York, 1986. t8
F T. Kroeh. “On Testing a Subroutine for the Numerical Integration Ordinary Dltferentlal equations,” J. ACM 20,1973, pp. 945-562.
19. J. A. Leavitt, “Methods and Applications vol. 20,1966, pp. 46-52. 20. R. E. Moore, 1966.
Interval
Analysis,
of Power Series,” Math. Comp.,
Prentice-Hall,
Englewood
21. R. E. Moore, Mathematical Elements of Scientific Rinehart and Winston, New York, 1975. 22. R. E. Moore, 1985.
Computational
Functional
23. L. B. Rail, Automatic Differentiation: Springer, New York, 1981.
of
Analysis, Techniques
Cliffs, N.J.,
Computing,
Wiley,
Holt,
New York,
and Applications,
24. L. F. Shampine, “Evaluation of a Test Set for Stiff ODE Solvers,” ACM Trans. Math. Soft., vol. 7, 1981, pp. 409-420. 25. R. E. Wengert, “A Simple Automafic Comm. ACM, vol. 9, 1964,443-464.
Derivative
Evaluation
Program”,