3. Mathematical Considerations

IDA solves the initial-value problem (IVP) for a DAE system of the general form

$$\prime \prime \prime$$ (3.1)

where y, y^$\prime$, and F are vectors in $\bf R^{N}_{}$, t is the independent variable, y^$\prime$ = dy/dt, and initial values y₀, y^$\prime$₀ are given. (Often t is time, but it certainly need not be.)

3.1 IVP solution

Prior to integrating a DAE initial-value problem, an important requirement is that the pair of vectors y₀ and y^$\prime$₀ are both initialized to satisfy the DAE residual F(t₀, y₀, y^$\prime$₀) = 0. For a class of problems that includes so-called semi-explicit index-one systems, IDA provides a routine that computes consistent initial conditions from a user's initial guess [4]. For this, the user must identify sub-vectors of y (not necessarily contiguous), denoted y_d and y_a, which are its differential and algebraic parts, respectively, such that F depends on y^$\prime$_d but not on any components of y^$\prime$_a. The assumption that the system is ``index one'' means that for a given t and y_d, the system F(t, y, y^$\prime$) = 0 defines y_a uniquely. In this case, a solver within IDA computes y_a and y^$\prime$_d at t = t₀, given y_d and an initial guess for y_a. A second available option with this solver also computes all of y(t₀) given y^$\prime$(t₀); this is intended mainly for quasi-steady-state problems, where y^$\prime$(t₀) = 0 is given. In both cases, IDA solves the system F(t₀, y₀, y^$\prime$₀) = 0 for the unknown components of y₀ and y^$\prime$₀, using Newton iteration augmented with a line search global strategy. In doing this, it makes use of the existing machinery that is to be used for solving the linear systems during the integration, in combination with certain tricks involving the step size (which is set artificially for this calculation). For problems that do not fall into either of these categories, the user is responsible for passing consistent values or risk failure in the numerical integration.

The integration method used in IDA is the variable-order, variable-coefficient BDF (Backward Differentiation Formula), in fixed-leading-coefficient form [1]. The method order ranges from 1 to 5, with the BDF of order q given by the multistep formula

$$\sum_{{i=0}}^{q} \alpha_{{n,i}}^{} \prime$$ (3.2)

where y_n and y^$\prime$_n are the computed approximations to y(t_n) and y^$\prime$(t_n), respectively, and the step size is h_n = t_n - t_n-1. The coefficients $\alpha_{{n,i}}^{}$ are uniquely determined by the order q, and the history of the step sizes. The application of the BDF (3.2) to the DAE system (3.1) results in a nonlinear algebraic system to be solved at each step:

$$\equiv \left(\vphantom{ t_n , y_n , h_n^{-1} \sum_{i=0}^q \alpha_{n,i}y_{n-i} }\right. \sum_{{i=0}}^{q} \alpha_{{n,i}}^{} \left.\vphantom{ t_n , y_n , h_n^{-1} \sum_{i=0}^q \alpha_{n,i}y_{n-i} }\right)$$ (3.3)

Regardless of the method options, the solution of the nonlinear system (3.3) is accomplished with some form of Newton iteration. This leads to a linear system for each Newton correction, of the form

J[y_n(m+1) - y_n(m)] = - G(y_n(m)) , (3.4)

where y_n(m) is the m-th approximation to y_n. Here J is some approximation to the system Jacobian

$${\frac{{\partial G}}{{\partial y}}} {\frac{{\partial F}}{{\partial y}}} \alpha {\frac{{\partial F}}{{\partial y^\prime}}}$$ (3.5)

where $\alpha$ = $\alpha_{{n,0}}^{}$/h_n. The scalar $\alpha$ changes whenever the step size or method order changes. The linear systems are solved by one of five methods:

• dense direct solver (serial version only),
• band direct solver (serial version only),
• diagonal approximate Jacobian solver,
• SPGMR = scaled preconditioned GMRES (Generalized Minimal Residual method) with restarts allowed,
• SPBCG = scaled preconditioned Bi-CGStab (Bi-Conjugate Gradient Stable method), or
• SPTFQMR = scaled preconditioned TFQMR (Transpose-Free Quasi-Minimal Residual method).

For the SPGMR, SPBCG, and SPTFQMR cases, preconditioning is allowed only on the left (see §3.2). Note that the direct linear solvers (dense and band) can only be used with serial vector representations.

In the process of controlling errors at various levels, IDA uses a weighted root-mean-square norm, denoted | ^. |_WRMS, for all error-like quantities. The multiplicative weights used are based on the current solution and on the relative and absolute tolerances input by the user, namely

W_i = 1/[RTOL ^. | y_i| + ATOL_i] . (3.6)

Because 1/W_i represents a tolerance in the component y_i, a vector whose norm is 1 is regarded as ``small.'' For brevity, we will usually drop the subscript WRMS on norms in what follows.

In the case of a direct linear solver (dense or banded), the nonlinear iteration (3.4) is a Modified Newton iteration, in that the Jacobian J is fixed (and usually out of date), with a coefficient $\bar{\alpha}$ in place of $\alpha$ in J. When using one of the Krylov methods SPGMR, SPBCG, or SPTFQMR as the linear solver, the iteration is an Inexact Newton iteration, using the current Jacobian (through matrix-free products Jv), in which the linear residual J$\Delta$y + G is nonzero but controlled. The Jacobian matrix J (direct cases) or preconditioner matrix P (SPGMR/SPBCG/SPTFQMR case) is updated when:

• starting the problem,
• the value $\bar{\alpha}$ at the last update is such that $\alpha$/$\bar{\alpha}$ < 3/5 or $\alpha$/$\bar{\alpha}$ > 5/3, or
• a non-fatal convergence failure occurred with an out-of-date J or P.

The above strategy balances the high cost of frequent matrix evaluations and preprocessing with the slow convergence due to infrequent updates. To reduce storage costs on an update, Jacobian information is always reevaluated from scratch.

The stopping test for the Newton iteration in IDA ensures that the iteration error y_n - y_n(m) is small relative to y itself. For this, we estimate the linear convergence rate at all iterations m > 1 as

$$\left(\vphantom{ \frac{\delta_m}{\delta_1} }\right. {\frac{{\delta_m}}{{\delta_1}}} \left.\vphantom{ \frac{\delta_m}{\delta_1} }\right)^{{\frac{1}{m-1}}}_{}$$

where the $\delta_{m}^{}$ = y_n(m) - y_n(m-1) is the correction at iteration m = 1, 2,.... The Newton iteration is halted if R > 0.9. The convergence test at the m-th iteration is then

$$\delta_{m}^{}$$ (3.7)

where S = R/(R - 1) whenever m > 1 and R$\le$0.9. The user has the option of changing the constant in the convergence test from its default value of 0.33. The quantity S is set to S = 20 initially and whenever J or P is updated, and it is reset to S = 100 on a step with $\alpha$ $\neq$ $\bar{\alpha}$. Note that at m = 1, the convergence test (3.7) uses an old value for S. Therefore, at the first Newton iteration, we make an additional test and stop the iteration if |$\delta_{1}^{}$| < 0.33 ^. 10^-4 (since such a $\delta_{1}^{}$ is probably just noise and therefore not appropriate for use in evaluating R). We allow only a small number (default value 4) of Newton iterations. If convergence fails with J or P current, we are forced to reduce the step size h_n, and we replace h_n by h_n/4. The integration is halted after a preset number (default value 10) of convergence failures. Both the maximum allowable Newton iterations and the maximum nonlinear convergence failures can be changed by the user from their default values.

When SPGMR, SPBCG, or SPTFQMR is used to solve the linear system, to minimize the effect of linear iteration errors on the nonlinear and local integration error controls, we require the preconditioned linear residual to be small relative to the allowed error in the Newton iteration, i.e., | P^-1(Jx + G)| < 0.05 ^. 0.33. The safety factor 0.05 can be changed by the user.

In the direct linear solver cases, the Jacobian J defined in (3.5) can be either supplied by the user or have IDA compute one internally by difference quotients. In the latter case, we use the approximation

$$\sigma_{j}^{} \prime \alpha \sigma_{j}^{} \prime \sigma_{j}^{}$$

$$\sigma_{j}^{} \sqrt{{U}} \left\{\vphantom{ \vert y_j\vert, \vert hy^\prime_j\vert,1/W_j }\right. \prime \left.\vphantom{ \vert y_j\vert, \vert hy^\prime_j\vert,1/W_j }\right\} \prime$$

where U is the unit roundoff, h is the current step size, and W_j is the error weight for the component y_j defined by (3.6). In the SPGMR/SPBCG/SPTFQMR case, if a routine for Jv is not supplied, such products are approximated by

$$\sigma \prime \alpha \sigma \prime \sigma$$

where the increment $\sigma$ is 1/| v|. As an option, the user can specify a constant factor that is inserted into this expression for $\sigma$.

During the course of integrating the system, IDA computes an estimate of the local truncation error, LTE, at the n-th time step, and requires this to satisfy the inequality

$$\leq$$

Asymptotically, LTE varies as h^q+1 at step size h and order q, as does the predictor-corrector difference $\Delta_{n}^{}$ $\equiv$ y_n - y_n(0). Thus there is a constant C such that

LTE = C$$\Delta_{n}^{}$$ + O(h^q+2) ,

and so the norm of LTE is estimated as | C| ^. |$\Delta_{n}^{}$|. In addition, IDA requires that the error in the associated polynomial interpolant over the current step be bounded by 1 in norm. The leading term of the norm of this error is bounded by $\bar{{C}}$|$\Delta_{n}^{}$| for another constant $\bar{{C}}$. Thus the local error test in IDA is

$$\bar{{C}} \Delta_{n}^{} \leq$$ (3.8)

A user option is available by which the algebraic components of the error vector are omitted from the test (3.8), if these have been so identified.

In IDA, the local error test is tightly coupled with the logic for selecting the step size and order. First, there is an initial phase that is treated specially; for the first few steps, the step size is doubled and the order raised (from its initial value of 1) on every step, until (a) the local error test (3.8) fails, (b) the order is reduced (by the rules given below), or (c) the order reaches 5 (the maximum). For step and order selection on the general step, IDA uses a different set of local error estimates, based on the asymptotic behavior of the local error in the case of fixed step sizes. At each of the orders q' equal to q, q - 1 (if q > 1), q - 2 (if q > 2), or q + 1 (if q < 5), there are constants C(q') such that the norm of the local truncation error at order q' satisfies

LTE(q') = C(q')|$$\phi$$(q' + 1)| + O(h^q'+2) ,

where $\phi$(k) is a modified divided difference of order k that is retained by IDA (and behaves asymptotically as h^k). Thus the local truncation errors are estimated as ELTE (q') = C(q')|$\phi$(q' + 1)| to select step sizes. But the choice of order in IDA is based on the requirement that the scaled derivative norms, | h^ky^(k)|, are monotonically decreasing with k, for k near q. These norms are again estimated using the $\phi$(k), and in fact

| h^q'+1y^(q'+1)| $$\approx$$ T(q') $$\equiv$$ (q' + 1)ELTE(q') .

The step/order selection begins with a test for monotonicity that is made even before the local error test is performed. Namely, the order is reset to q' = q - 1 if (a) q = 2 and T(1) $\leq$ T(2)/2, or (b) q > 2 and max{T(q - 1), T(q - 2)} $\leq$ T(q); otherwise q' = q. Next the local error test (3.8) is performed, and if it fails, the step is redone at order q $\leftarrow$ q' and a new step size h'. The latter is based on the h^q+1 asymptotic behavior of ELTE(q), and, with safety factors, is given by

$$\eta$$ = h'/h = 0.9/[2 ELTE(q)]^1/(q+1) .

The value of $\eta$ is adjusted so that 0.25 $\leq$ $\eta$ $\leq$ 0.9 before setting h $\leftarrow$ h' = $\eta$h. If the local error test fails a second time, IDA uses $\eta$ = 0.25, and on the third and subsequent failures it uses q = 1 and $\eta$ = 0.25. After 10 failures, IDA returns with a give-up message.

As soon as the local error test has passed, the step and order for the next step may be adjusted. No such change is made if q' = q - 1 from the prior test, if q = 5, or if q was increased on the previous step. Otherwise, if the last q + 1 steps were taken at a constant order q < 5 and a constant step size, IDA considers raising the order to q + 1. The logic is as follows: (a) If q = 1, then reset q = 2 if T(2) < T(1)/2. (b) If q > 1 then

• reset q $\leftarrow$ q - 1 if T(q - 1) $\leq$ min{T(q), T(q + 1)};
• else reset q $\leftarrow$ q + 1 if T(q + 1) < T(q);
• leave q unchanged otherwise [then T(q - 1) > T(q) $\leq$ T(q + 1)].

In any case, the new step size h' is set much as before:

$$\eta$$ = h'/h = 1/[2 ELTE(q)]^1/(q+1) .

The value of $\eta$ is adjusted such that (a) if $\eta$ > 2, $\eta$ is reset to 2; (b) if $\eta$ $\leq$ 1, $\eta$ is restricted to 0.5 $\leq$ $\eta$ $\leq$ 0.9; and (c) if 1 < $\eta$ < 2 we use $\eta$ = 1. Finally h is reset to h' = $\eta$h. Thus we do not increase the step size unless it can be doubled. See [1] for details.

IDA permits the user to impose optional inequality constraints on individual components of the solution vector y. Any of the following four constraints can be imposed: y_i > 0, y_i < 0, y_i $\geq$ 0, or y_i $\leq$ 0. The constraint satisfaction is tested after a successful nonlinear system solution. If any constraint fails, we declare a convergence failure of the Newton iteration and reduce the step size. Rather than cutting the step size by some arbitrary factor, IDA estimates a new step size h' using a linear approximation of the components in y that failed the constraint test (including a safety factor of 0.9 to cover the strict inequality case). These additional constraints are also imposed during the calculation of consistent initial conditions.

Normally, IDA takes steps until a user-defined output value t = t_out is overtaken, and then computes y(t_out) by interpolation. However, a ``one step'' mode option is available, where control returns to the calling program after each step. There are also options to force IDA not to integrate past a given stopping point t = t_stop.

3.2 Preconditioning

When using a Newton method to solve the nonlinear system (3.4), IDA makes repeated use of a linear solver to solve linear systems of the form J$\Delta$y = - G. If this linear system solve is done with one of the scaled preconditioned iterative linear solvers, these solvers are rarely successful if used without preconditioning; it is generally necessary to precondition the system in order to obtain acceptable efficiency. A system Ax = b can be preconditioned on the left, on the right, or on both sides. The Krylov method is then applied to a system with the matrix P^-1A, or AP^-1, or P_L^-1AP_R^-1, instead of A. However, within IDA, preconditioning is allowed only on the left, so that the iterative method is applied to systems (P^-1J)$\Delta$y = - P^-1G. Left preconditioning is required to make the norm of the linear residual in the Newton iteration meaningful; in general, | J$\Delta$y + G| is meaningless, since the weights used in the WRMS-norm correspond to y.

In order to improve the convergence of the Krylov iteration, the preconditioner matrix P should in some sense approximate the system matrix A. Yet at the same time, in order to be cost-effective, the matrix P should be reasonably efficient to evaluate and solve. Finding a good point in this tradeoff between rapid convergence and low cost can be very difficult. Good choices are often problem-dependent (for example, see [2] for an extensive study of preconditioners for reaction-transport systems).

Typical preconditioners used with IDA are based on approximations to the Newton iteration matrix of the systems involved; in other words, P $\approx$ ${\frac{{\partial F}}{{\partial y}}}$ + $\alpha$${\frac{{\partial F}}{{\partial y^\prime}}}$, where $\alpha$ is a scalar inverse proportional to the integration step size h. Because the Krylov iteration occurs within a Newton iteration and further also within a time integration, and since each of these iterations has its own test for convergence, the preconditioner may use a very crude approximation, as long as it captures the dominant numerical feature(s) of the system. We have found that the combination of a preconditioner with the Newton-Krylov iteration, using even a fairly poor approximation to the Jacobian, can be surprisingly superior to using the same matrix without Krylov acceleration (i.e., a modified Newton iteration), as well as to using the Newton-Krylov method with no preconditioning.

3.3 Rootfinding

The IDA solver has been augmented to include a rootfinding feature. This means that, while integrating the Initial Value Problem (3.1), IDA can also find the roots of a set of user-defined functions g_i(t, y, y') that depend on t, the solution vector y = y(t), and its t -derivative y'(t). The number of these root functions is arbitrary, and if more than one g_i is found to have a root in any given interval, the various root locations are found and reported in the order that they occur on the t axis, in the direction of integration.

Generally, this rootfinding feature finds only roots of odd multiplicity, corresponding to changes in sign of g_i(t, y(t), y'(t)), denoted g_i(t) for short. If a user root function has a root of even multiplicity (no sign change), it will probably be missed by IDA. If such a root is desired, the user should reformulate the root function so that it changes sign at the desired root.

The basic scheme used is to check for sign changes of any g_i(t) over each time step taken, and then (when a sign change is found) to home in on the root (or roots) with a modified secant method [11]. In addition, each time g is computed, IDA checks to see if g_i(t) = 0 exactly, and if so it reports this as a root. However, if an exact zero of any g_i is found at a point t, IDA computes g at t + $\delta$ for a small increment $\delta$, slightly further in the direction of integration, and if any g_i(t + $\delta$) = 0 also, IDA stops and reports an error. This way, each time IDA takes a time step, it is guaranteed that the values of all g_i are nonzero at some past value of t, beyond which a search for roots is to be done.

At any given time in the course of the time-stepping, after suitable checking and adjusting has been done, IDA has an interval (t_lo, t_hi] in which roots of the g_i(t) are to be sought, such that t_hi is further ahead in the direction of integration, and all g_i(t_lo) $\neq$ 0. The endpoint t_hi is either t_n, the end of the time step last taken, or the next requested output time t_out if this comes sooner. The endpoint t_lo is either t_n-1, or the last output time t_out (if this occurred within the last step), or the last root location (if a root was just located within this step), possibly adjusted slightly toward t_n if an exact zero was found. The algorithm checks g at t_hi for zeros and for sign changes in (t_lo, t_hi). If no sign changes are found, then either a root is reported (if some g_i(t_hi) = 0) or we proceed to the next time interval (starting at t_hi). If one or more sign changes were found, then a loop is entered to locate the root to within a rather tight tolerance, given by

$$\tau$$ = 100*U*(| t_n| + | h|)   (U = unit roundoff) .

Whenever sign changes are seen in two or more root functions, the one deemed most likely to have its root occur first is the one with the largest value of | g_i(t_hi)|/| g_i(t_hi) - g_i(t_lo)|, corresponding to the closest to t_lo of the secant method values. At each pass through the loop, a new value t_mid is set, strictly within the search interval, and the values of g_i(t_mid) are checked. Then either t_lo or t_hi is reset to t_mid according to which subinterval is found to have the sign change. If there is none in (t_lo, t_mid) but some g_i(t_mid) = 0, then that root is reported. The loop continues until | t_hi - t_lo| < $\tau$, and then the reported root location is t_hi.

In the loop to locate the root of g_i(t), the formula for t_mid is

t_mid = t_hi - (t_hi - t_lo)g_i(t_hi)/[g_i(t_hi) - $$\alpha$$g_i(t_lo)] ,

where $\alpha$ a weight parameter. On the first two passes through the loop, $\alpha$ is set to 1, making t_mid the secant method value. Thereafter, $\alpha$ is reset according to the side of the subinterval (low vs high, i.e. toward t_lo vs toward t_hi) in which the sign change was found in the previous two passes. If the two sides were opposite, $\alpha$ is set to 1. If the two sides were the same, $\alpha$ is halved (if on the low side) or doubled (if on the high side). The value of t_mid is closer to t_lo when $\alpha$ < 1 and closer to t_hi when $\alpha$ > 1. If the above value of t_mid is within $\tau$/2 of t_lo or t_hi, it is adjusted inward, such that its fractional distance from the endpoint (relative to the interval size) is between .1 and .5 (.5 being the midpoint), and the actual distance from the endpoint is at least $\tau$/2.