On Implicit Taylor Series Methods for Stiff ODEs

G. Kirlinger  G. F. Corliss^{∗}
Institut für Angewandte  Mathematics and Computer Science Division
und Numerische Mathematik Argonne National Laboratory
Technische Universität Wien  9700 S. Cass Avenue
Wiedner Hauptstrasse 8–10  Argonne, IL 60439–4801 USA

A-1040 Wien, AUSTRIA

Abstract

Several versions of implicit Taylor series methods (ITSM) are presented and evaluated. Criteria for the approximate solution of ODEs via ITSM are given. Some ideas, motivations, and remarks on the inclusion of the solution of stiff ODEs are outlined.

1 Introduction

One approach for the validated solution of initial value problems for ODEs is built upon a Taylor series method. Moore [20, 21], Rall [22], Adams [1], Lohner [17], and Eijgenraam [14] were the first to use Taylor series methods to enclose the solutions of ODEs. One of the main advantages of using a Taylor series method is the simple representation of the local discretization error. For a Taylor series method, the discretization error is the remainder term of the series, which can easily be bounded by using automatic differentiation to generate as many derivatives as required. Further, the order of the Taylor series can also easily be adjusted to the needed accuracy.

∗Supported in part by the National Science Foundation grant No. CCR–8802429, by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S.

Department of Energy, under Contract W-31-109-Eng-38, and through NSF Cooperative Agreement No. CCR-8809615. On leave from Marquette University, Milwaukee, WI 53233 USA.

This paper emphasizes approximation methods based on Taylor series, but we also give some initial ideas, motivations, and remarks on the inclusion of the solution of stiff ODEs (Section 5). We compare several numerical methods based on Taylor series to evaluate their suitability for use as methods for the verified inclusion of stiff differential equations.

This is a report on work in progress. Some numerical experiments have been done, but much more experimental and theoretical work remains to be done. Another manuscript [7] is in preparation describing in detail the relationship of the implicit Taylor series methods described here to implicit methods based on Padé approximations.

2 Stiffness

In this paper, we present initial steps towards algorithms for solving stiff systems of ODEs based on Taylor series. We begin with a brief survey of theory and point algorithms for stiff problems. The modeling of evolution, or growth, processes in many applied sciences often leads to stiff ordinary differential equations. The initial value problem

y′ = λy, y(0) = y_{0} with λ ≪ 0  (2.1)

is a simple example of a stiff problem. It shows the following characteristic behavior:

(i) The Lipschitz constant L = |λ| of the right-hand side is very large.

(ii) For y_{0} = 0, the solution y(t) ≡ 0 is smooth.

(iii) For y_{0} ≠ 0, the solution y(t) = e^{λt}y_{0} is rapidly decaying and very unsmooth in the sense that derivatives are large during the initial “transient phase.” Away from t = 0, the solution becomes smooth very quickly and tends towards the smooth solution (“smooth phase”).

In most practical situations, rapidly decaying components (corresponding to λ ≪ 0) occur together with smooth (nonstiff) components. Such a system quickly tends towards an “equilibrium,” that is, to a smooth solution. Rapid variations occur only while the state of equilibrium has not yet been reached or when the system is switched from one state to another (e.g., by nonlinear effects).

Although it is difficult to give a mathematically rigorous definition of stiffness, we call a system

y′ = f(y)  (2.2)
y(0) = y_{0}

stiff if its Jacobian f_{y} (in a neighborhood of the solution) has eigenvalues λ_{i} with Re(λ_{i}) ≪ 0, in addition to eigenvalues of moderate size. (The autonomous notation of (2.2) is only for technical simplification of the formulas in sections 3 to 5.)

Stiff systems are considered difficult because explicit numerical methods designed for nonstiff problems are forced to use very small steps. If “normal” steps are used, then perturbations in the computed solution are amplified by the influence of the Lipschitz constant L ≫ 0. In order to retain the stability of the true solution in the computed solution, the step must be very small. It is not possible to use a step size that is adjusted to the smoothness of the solution sought. That is why authors as early as 1928 [8] or 1947 [9] were led to consider implicit methods in which the approximate solution y_{i+1} at t = t_{i+1} is given by the solution to some nonlinear system. Many implicit methods allow step sizes appropriate to the smoothness of the solution.

Historically, the first theoretical concept especially suited for the assessment of numerical methods for stiff problems was A-stability [11]. A-stability means that computed numerical approximations to decreasing components are also decreasing. This analysis is based on the model linear constant coefficient problem y′ = Ay. The concept of B-stability [15] is the basis for a general convergence theory for nonlinear stiff problems. Nevertheless, we use the linear concept A-stability as a first criterion for the assessment of implicit procedures based on Taylor series methods.

The general algorithm for implicit methods for stiff systems is as follows:

Initialize Loop for each integration step

Guess step size

Solve some nonlinear system for y_{i+1} Estimate error

Accept or reject step

Well-known examples of methods of this class that are suitable for the integration of stiff systems are backward differentiation formulas (BDFs) [10] and implicit Runge-Kutta methods [4] and [12]. Our contribution is to use Taylor series to formulate the nonlinear system for y_{i+1}.

Any numerical one-step method applied to (2.1) reduces to

y_{i+1} = R(z)y_{i}  z = hλ,  (2.3)

where R : C → C is a polynomial or (in the case of implicit methods) a rational function with real coefficients. R(z) is called the stability function of the method.

Definition: A numerical method applied to (2.1) as well as its stability function R(z) is called A-stable if the left half-plane {z : Re(z) ≤ 0} is contained in the region {z : |R(z)| ≤ 1}.

On the one hand, the definition of A-stability is too weak: only linear problems with constant coefficients are covered. On the other hand, the definition is too strong: many methods that perform well in practice are not A-stable. Thus we are led to the following generalization.

Definition: A numerical method is A(α)-stable if the sector {z : |arg(−z)| ≤ α, z ≠ 0} is contained in the stability region {z : |R(z)| ≤ 1}.

In contrast to A-stable methods, there exists little nonlinear convergence theory in the literature for methods that are only A(α)-stable. A step in this direction for BDFs was given by Lubich [19].

We were led to consider implicit Taylor series methods for the validated solution of initial value problems for stiff ODEs because explicit Taylor series methods that have proven effective for nonstiff systems [6, 2] cannot be expected to handle highly stiff systems successfully and because Lohner’s program AWA [18]) has proven so successful using Taylor series methods for computing inclusions.

3 Implicit Taylor Series Methods (ITSM)

The first formulations of ITSM are due to Chang and Corliss [5] and to Stetter [25]. Let us start by recalling explicit Taylor series methods. Let y_{i} be an approximation for the solution of (2.2) at t = t_{i}. The explicit Taylor series method uses recurrence relations derived from the ODE to generate the series for y(t) expanded at t_{i}. Then

y_{i+1} := y_{i} + ∑_{j=1}^{k} f_{j−1}(y_{i}) \frac{h^{j}}{j!} ,  (3.4)

where f_{j−1}(y) := \frac{d^{j}y}{dt^{j}} for j = 0, …, k − 1, and h := t_{i+1} − t_{i} is taken as large as possible consistent with error control constraints.

The simplest form of the ITSM is as follows:

Initialize Loop for each integration step

Guess step size

Loop Guess or improve y_{i+1} Generate series for y expanded at t = t_{i+1} Exit loop when y_{i} = y_{i+1} + ∑_{j=1}^{k} f_{j−1}(y_{i+1}) \frac{(t_{i}−t_{i+1})^{j}}{j!} is satisfied within tolerance

Estimate error

Accept or reject step

The series at t = t_{i+1} is generated by using exactly the same recurrence relations as in the explicit method. That is, each f_{j−1}(y_{i+1}) is some nonlinear function of y_{i+1}. The difference is that the “initial condition” y_{i+1} is determined by a Newton iteration, rather than being given by analytic continuation at the previous step. See Figure 1.

The Newton iteration for the equation

y_{i} = y_{i+1} + ∑_{j=1}^{k} f_{j−1}(y_{i+1}) \frac{(t_{i} − t_{i+1})^{j}}{j!}  (3.5)

Figure 1: Implicit Taylor series method

requires

\frac{∂}{∂y_{i+1}} (y_{i+1} + ∑_{j=1}^{k} f_{j−1}(y_{i+1}) \frac{(t_{i} − t_{i+1})^{j}}{j!} ) = 1I+∑_{j=1}^{k} \frac{∂}{∂y_{i+1}} (f_{j−1}(y_{i+1})) \frac{(t_{i} − t_{i+1})^{j}}{j!} ,

(3.6)

which is computed from the solution of the variational equation U′ = \frac{∂f}{∂y} U [17].

The stability function R(z) corresponding to this method is a (0, k)−Padé approximation to the exponential exp(z) [12], where k is the degree of the Taylor series. It is well known that the (m, k)−Padé approximation is the unique rational function with numerator and denominator of degree m and k, respectively, which approximates exp(z) to O(z^{k+m+1}) as z → 0.

An (m, k)−Padé approximation to exp(z) is A-stable if and only if k − 2 ≤ m ≤ k. That is, the diagonal and two subdiagonals in the Padé scheme are A-stable (see [12] and [13]).

In our case, only implicit Taylor series of degree 1 and 2 lead to A-stable ITSM. For degree k = 1, the resulting method is the implicit Euler method, whose stability and convergence properties are well known. For degree k = 2, the resulting method corresponds to a special Runge-Kutta scheme, Lobatto IIIC with the number of stages s = 2. The corresponding stability function is a (s − 2, s)−Padé approximation. Schneid [24] showed that the Lobatto Figure 2: Stability regions of the (0, k)−Padé approximation for k =1, 2, 3, 4, 5, 6

IIIC method for s = 2 is not only A-stable, but also B-convergent of order 2 under some reasonable step-size restrictions.

The corresponding stability regions for the (0, k)−Padé approximation (see Figure 2) were drawn by using the software product S [3]. Sand and Østerby [23] give stability regions for certain Runge-Kutta methods that are qualitatively similar to Figure 2.

4 Implicit σ Taylor Series Methods (IσTSM)

The lack of A-stability of the simple ITSM for higher orders suggests that we try to increase the degree of the numerator in the corresponding rational stability function.

The ITSM matches the previously computed y_{i}, with the value obtained by expanding the solution at t_{i+1}. The IσTSM generates the Taylor series at both t_{i} and t_{i+1}. The condition that the two series agree at t_{i} + σ(t_{i+1} − t_{i}), σ ∈ (0, 1), provides the nonlinear equation for y_{i+1} (see Figure 3).

Figure 3: Implicit σ Taylor series method

The test

y_{i} ≈ y_{i+1} + ∑_{j=1}^{k} f_{j−1}(y_{i+1}) \frac{(t_{i} − t_{i+1})^{j}}{j!}  (4.7)

in the ITSM is now replaced by the test

y_{i} + ∑_{j=1}^{m} f_{j−1}(y_{i}) \frac{σ^{j}(t_{i+1} − t_{i})^{j}}{j!} ≈ y_{i−1} + ∑_{j=1}^{k} f_{j−1}(y_{i+1}) \frac{(σ − 1)^{j}(t_{i+1} − t_{i})^{j}}{j!} .

(4.8)

In general, m need not be equal to k. The IσTSM requires almost no work beyond that required by the ITSM. The series at t_{i+1} must be recomputed for each new iterate y_{i+1}, as in the ITSM described in Section 3. The series at t_{i}, however, is computed only once per step, and this computation has already been done at the end of the previous step. Here, σ is a tuning parameter of the method. The case σ = 1 represents an explicit Taylor series, whereas σ = 0 leads to the fully implicit form of Section 3. The case σ = \frac{1}{2} is the unique choice of σ for which the resulting method is A-stable, for any order m = k. Unfortunately, the resulting stability function is not a Padé approximation if the orders m, k of the Taylor series are higher than 1. Hence, the maximal achievable order of the local truncation error is reduced compared to the maximal possible order for the Padé approximation.

Remark: The IσTSM with k = m = 1 and σ = \frac{1}{2} is the well-known implicit trapezoidal rule, which is A-stable, but not B-stable.

In a forthcoming paper [7], a class of high-order stiff ordinary differential equations solvers based on Padé approximations is introduced and analyzed. In this approach, the advantages of Taylor series methods mentioned above are maintained and combined with highest possible order of the local truncation error. A general nonlinear stability analysis remains to be done.

5 Implications for Interval Methods

The first task when developing an algorithm for validated inclusions of solutions to stiff problems is to understand an interval analog of the implicit Euler scheme y_{i+1} = y_{i} + hf(y_{i+1}). Let ỹ_{i} represent an approximate solution, and let [e_{i}] be an interval inclusion for the corresponding error. We represent an interval-valued solution at t_{i} as [y_{i}] = ỹ_{i} + [e_{i}]. The implicit Euler method becomes

e_{n+1} − hf_{y}(ỹ_{n+1} + [0, 1] e_{n+1}) · e_{n+1} = e_{n} − \frac{h^{2}}{2} f_{y}(ȳ) · f(ȳ),  (5.9)

for e_{n} ∈ [e_{n}], and ȳ contained in an a-priori inclusion Ȳ_{n+1}; the approximate solution ỹ_{n+1} at t_{n+1} is chosen such that ỹ_{n+1} = ỹ_{n} + hf(ỹ_{n+1}). This is a nonlinear system of interval equations. In the special case of a linear constant coefficient problem y′ = Ay, the corresponding system reduces to a linear system with interval right-hand side. The methods presented in this paper can be viewed as higher-order generalizations of this simplest implicit scheme.

Another difficulty in the case of stiffness is the inevitably large Lipschitz constant L of the right-hand f(y). This would cause a severe and unrealistic step-size restriction whenever the Picard-Lindelöf existence theorem is applied to get an a-priori inclusion. Subtle algorithms have to be developed.

During the preparation of this manuscript, we received Kreuser’s diploma thesis [16], which presents an alternative approach to computing inclusions of solutions for stiff ODEs. Kreuser’s approach transforms the original stiff system (2.2) into the new system

y′ = Ay + (f(y) − Ay),

where A is a local approximation of the Jacobian f_{y}(y). The a-priori inclusion is obtained by symbolic computation of the matrix exponential function, thus avoiding the severe step-size restriction inherent in Lohner’s explicit approach [17]. To get an inclusion [y_{i}] = ỹ_{i} + [e_{i}] for the solution at t_{i}, the approximate solution ỹ_{i} is obtained by integrating the linear ODE and applying a Taylor series method to the nonlinear part. The inclusion [e_{i}] is obtained by iterating an interval integral operator. The method of Kreuser successfully handles stiff problems of the special form y′ = Ay + g(t, y), where g(t, y) has a small Lipschitz constant w.r.t. y. This excludes strongly nonlinear models.

Acknowledgments

The authors thank Prof. H. J. Stetter and W. Auzinger for many valuable discussions. Special thanks are also due to an anonymous gardener from the Viennese Stadtpark who found and returned the first author’s backpack after it had been stolen; it contained all written records and source material for the work reported here.

References

[1] E. Adams, Periodic solutions: Enclosure, verification, application, in Computer Arithmetic and Self Validating Numerical Methods, Proceedings of SCAN 1989 (Basel), Academic Press, Boston, 1990, 199–245.

[2] D. Barton, I. M. Willers, R. V. M. Zahar, The automatic solution of ordinary differential equations by the method of Taylor series, Comp. J., 14 (1971), 243–248.

[3] R. A. Becker, J. M. Chambers, A. R. Wilks, The New S Language: A Programming Environment for Data Analysis and Graphing, Wadsworth and Brooks/Cole, Pacific Grove, Calif., 1988.

[4] J. C. Butcher, Implicit Runge-Kutta processes, Math. Comp., 18 (1964) 50–64.

[5] Y. F. Chang, G. F. Corliss, Personal communication.

[6] G. F. Corliss, Y. F. Chang, Solving ordinary differential equations using Taylor series, ACM Trans. Math. Software, 8 (1982), 114–144.

[7] G. F. Corliss, A. Griewank, F. Potra, G. Kirlinger, H. J. Stetter, Personal communication.

[8] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann., 100 (1928), 32–74.

[9] J. Crank, P. Nicholson, A practical method for numerical integration of solutions of partial differential equations of heat conduction type, Proc. Cambridge Philos. Soc., 43 (1947), 50 ff.

[10] C. F. Curtis, J. O. Hirschfelder, Integration of stiff equations, Proc. Nat. Acad. Sci., 38 (1952), 235–243.

[11] G. Dahlquist, A special stability problem for linear multistep methods, BIT, 3 (1963), 27–43.

[12] K. Dekker, J. G. Verwer, Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984.

[13] B. L. Ehle, A-stable methods and Padé approximations to the exponential, SIAM J. Math. Anal. 4 (1973), 671–680.

[14] P. Eijgenraam, The Solution of Initial Value Problems Using Interval Arithmetic, Stichting mathematisch Centrum Amsterdam, Mathematical Center Tracts No. 144, 1981.

[15] R. Frank, J. Schneid, Ch. W. Überhuber, The concept of B-convergence, SIAM J. Numer. Anal., 18 (1981), 753–780.

[16] A. Kreuser, Zur Lösungseinschließung bei steifen Differentialgleichungen durch eine Verbindung von Computeralgebra und Intervallarithmetik, Diplomarbeit (Universität Bonn), 1991.

[17] R. Lohner, Einschließung der Lösung gewöhnlicher Anfangs- und Randwertaufgaben und Anwendungen, Dissertation (Universität Karlsruhe), 1988.

[18] R. Lohner, AWA, Software product in FORTRAN-SC for the inclusion of the solution of ODEs, Karlsruhe, 1989.

[19] Ch. Lubich, On the convergence of multistep methods for nonlinear stiff differential equations, Numer. Math. 58 (1991) 839–853.

[20] R. E. Moore, Interval Analysis, Prentice-Hall, Englewood Cliffs, N.J., 1966.

[21] R. E. Moore, Methods and Applications of Interval Analysis, SIAM, Philadelphia, 1979.

[22] L. B. Rall, Automatic Differentiation: Techniques and Applications, Lecture Notes in Comp. Sci. No. 120, Springer-Verlag, Berlin, 1981.

[23] J. Sand, O. Østerby, Regions of absolute stability, DAIMI PB-102, Computer Science Department, Aarhuis University, September 1979.

[24] J. Schneid, B-convergence of Lobatto IIIC formulas, Num. Math., 51 (1987), 229–235.

[25] H. J. Stetter, Validated solution of initial value problems for ODE, in Computer Arithmetic and Self Validating Numerical Methods, Proceedings of SCAN 1989 (Basel), Academic Press, Boston, 1990, 171–187.