Consistent Initial Condition Calculation for Differential-Algebraic Systems.
In this paper we describe a new algorithm for the calculation of consistent initial conditions for a class of systems of differential-algebraic equations which includes semi-explicit index-one systems.  We consider initial condition problems of two types---one where the differential variables are specified, and one where the derivative vector is specified.  The algorithm requires a minimum of additional information from the user.  We outline the implementation in a general-purpose solver DASPK for differential-algebraic equations, and present some numerical experiments which illustrate its effectiveness.
Introduction
. This paper is concerned with the calculation of initial conditions for systems
of differential-algebraic equations (DAEs). We write the DAE system in the general form
G(t;
G, y, and y 0 are N-dimensional vectors. The initial value problem for this system is the
problem of finding a solution that satisfies a consistent set of initial conditions y(t 0
. Two software packages have been written for solving initial value problems for the DAE system
(1.1)-DASSL [2], and an extension of it called DASPK [7]. Both use variable-order, variable-
stepsize backward differentiation formulas. DASSL solves the linear system at each time step by
dense or banded direct linear system methods. In DASPK, the linear systems that arise at each
time step are solved with either direct linear system methods, or with a preconditioned Krylov
iterative method, namely GMRES [16]. For large-scale systems, the iterative method combined
with a suitable preconditioner can be quite effective.
When using either of the solvers DASSL or DASPK, the integration must be started with a
consistent set of initial conditions y 0 and y 0
. Consistency requires, in particular, that G(t 0 ; y
Usually, not all of the components of y 0 and y 0
are known directly from the original problem
specification. The problem of finding consistent initial values can be a challenging task. The present
DASSL and DASPK solvers offer an option for finding consistent y 0
0 from a given initial y 0 , by taking
a small artificial step with the Backward Euler method. However, initialization problems do not
always arise in this form, and even for the intended problem type, that technique is not always
successful. In any case it is unsatisfactory in that it produces values at
than at . In this paper, we propose an alternative procedure for a class of DAE problems.
We will show that this method, in combination with the modified Newton methods of DASSL or
the Newton-Krylov methods of [5] and [6], yields an algorithm which converges nearly as rapidly
as the underlying Newton or Newton-Krylov method. The new method is very convenient for the
user, because it makes use of the Jacobian or preconditioner matrices which are already required in
DASSL or DASPK. The consistent initialization problem has been studied in [1, 8, 9, 12, 13, 15, 14]
for more general DAE systems than we will consider here. However, none of these algorithms offer
the simplicity or convenience of the approach which we develop here.
This research was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy
Research, U.S. Dept. of Energy, by Lawrence Livermore National Laboratory under contract W-7405-ENG-48.
y Center for Computational Sciences & Engineering, L-316, Lawrence Livermore National Laboratory, Livermore,
California 94551.
z Department of Computer Science, University of Minnesota, Minneapolis, MN 55455. The work of this author
was partially supported by ARO contract number DAAL03-92-G-0247, DOE contract number DE-FG02-92ER25130,
NIST contract number 60NANB2D1272, and by the Minnesota Supercomputer Institute.
The class of problems that we consider is a generalization of semi-explicit index-one DAE
systems. Semi-explicit index-one DAE systems are characterized as follows. The dependent variable
vector y can be split into a vector u of size N d , called the differential variables, and a vector v of
size N a , called the algebraic variables, such that the equations have the form
in which g v = @g=@v is a nonsingular square matrix. We will be concerned with the initialization
problem of finding the initial value v 0 of v when the initial value u 0 for u is specified.
We can easily generalize the class of problems considered, to those where the ODE subsystem
for u may be implicit. Thus we consider systems of the form
where with the matrix f u being square and nonsingular.
We will continue to refer to problems of the form (1.3), with f u 0 and g v nonsingular, as semi-explicit
index-one, even though they may be not be explicit in u 0 . In fact, our main initialization technique
applies to an even more general class of problems, as we will explain later.
We also consider a second type of initialization problem, in which the initial derivatives are
specified but all of the dependent variables are unknown. That is, we must solve for y 0 given y 0
For example, beginning the DAE solution at a steady state corresponds to specifying y 0
problem does not involve a split of the y vector into differential and algebraic parts or a semi-explicit
form for the equations.
In later sections, we will refer to these two problems as Initialization Problem 1 and Initialization
Problem 2.
2. The Basic Method. The central idea of this paper is to solve both of these initial condition
problems with the help of mechanisms already in place for the solution of the DAE system itself,
rather than requiring the user to perform a special computation for it.
Consider first Initialization Problem 1 for the semi-explicit index-one system (1.3), where v
is to be determined, given u at the initial point We expand this problem to
include the calculation of u 0
Thus we can form a nonlinear system in the N-vector
namely
This general approach of solving the expanded problem has also been used in [12]. A Newton
iteration for the solution of F would require the Jacobian matrix
f
By assumption, this matrix is nonsingular, at least in a neighborhood of the desired solution.
In the course of integrating a DAE system with DASSL or DASPK, the user must call upon
one of several linear system algorithms to solve N \Theta N linear systems at every time step. These
arise from a Newton-like method for solving algebraic systems G(t; y, where a is
a vector containing past values, and c is a constant, set by the solver, that is inversely proportional
to the stepsize h. Thus the linear systems have the form
in which R is a residual vector, \Deltay is a correction to y, and the matrix J is the DAE system
iteration matrix
@y
The user is encouraged to supply an approximation to J , for use either as the Newton matrix itself
in the case of direct methods, or as a preconditioner in the case of a Krylov method. In the direct
case, J is generated by difference quotient approximations if not supplied by the user. In the case
of the system (1.3), we have
f
f
cf
In order to make use of J in solving F we pick an artificial stepsize h, and set
in (2.5). Then, to recover the block f u 0 , we rescale the first block-column of J by h, using the
scaling matrix
I a
where I d and I a are the identity matrices of size N d and N a , respectively. Thus we consider the
matrix
hg
evaluated at
. Note that -
f
Thus if h is small in some appropriate sense, we can expect that -
J will be a good approximation
to F 0 (x).
The proposed initialization procedure is to carry out a Newton-like iteration with corrections
Each iteration will call on the linear system solution procedure that is to be used later in solving
the DAE system itself. It will also require information about which components of y are differential
and which are algebraic, in order to apply the correction \Deltax to the vectors y and y 0 . But otherwise,
the procedure requires no additional information or methodology. Upon convergence, we have all
components of y(t 0 ), and we have the components of y corresponding to u 0
0 , the derivatives of
the differential variables. The remaining components of y corresponding to v 0
simply be
set to zero, as the integration procedure is insensitive to these (since v 0 does not appear in (1.3)) 1 ,
1 Although the BDF formulas do not depend on these values, the error test at the end of the first step depends
on them, unless the v variables are excluded from the error test (which is an option in the code).
and the first time step will produce accurate values for them. The next two sections will present
this procedure in a more formal manner and prove convergence for it.
For Initialization Problem 2, we are given the initial value of y 0 and must compute the initial
y. In this case, we are simply interested in solving for in the system
with y 0
We assume that this problem is well-posed as provided by the user, so that F 0
G y is nonsingular in a neighborhood of the solution, including the initial guess supplied. As in the
first problem, we will call for the user to supply the DAE iteration matrix J , but this time we set
so that the matrix involved is simply there is no stepsize h. We then proceed
with Newton iterations using J , with corrections
\Deltay
Finally, we remark that an extension of these ideas to Hessenberg index-2 DAE systems is
under way by the authors.
3. More General Problems. In the case of Initialization Problem 1, the full scope of problems
for which the above idea can be applied is more general than systems of the form (1.3).
0 and satisfies F The true
Jacobian is F 0 and is approximated by -
suitably
small.
3.2. Permuted variables. We wish to generalize the system (1.3) first by dropping the requirement
that the differential and algebraic components are separated into blocks in y. Thus we
assume there is a permutation matrix P of size N such that
while in terms of z the system function G has the form
The vector of unknowns in terms of z is
and in terms of y it is
In order to express the objective function F (x) precisely, we introduce two projections in R N ,
I d 0
I a
These satisfy the help of these projections, we can write
I a
I a
and
P a y
I a
That is, P a y a x is a vector containing the unknown algebraic components (corresponding to
Similarly,
vector containing the unknown derivative components (corresponding to u 0
where. The known components of y 0 comprise a vector
and so we can write y x. The fact that G does not depend on v 0 means that G(t;
Thus the system to be solved in x is
In the same way, we need not require that the components of G are blocked as in (3.2). Thus
we allow a permutation Q in the components of G, such that QG has that blocked form. Then
of course so does QF . However, we can work with G (hence F ) in its original ordering, in the
initialization procedure.
3.3. Implicit constraints. To generalize further the form of the problem we can solve, suppose
that, after permuting the y vector to z = Py, the DAE system function G has the form
G(t; having full rank N d :
This class of problems generalizes (1.3) and (3.2) in that the algebraic constraints, g(t; u;
(1.3), need not be identified explicitly.
For systems of this form, we can again define the projection matrices P d and P a by (3.5), and
it again is true that
G(t;
We define the vectors w and x by (3.3) and (3.4) as before, and the system to be solved is now
with y before. The Jacobian of this function is
Knowing that H u 0 has full rank, there is a nonsingular matrix M such that
MH
A 1'
with A 1 a square nonsingular matrix of size N d . Then
and the fact that the index is one implies that B 2 is nonsingular near the consistent initial value
(see x2.5 of [2]). Thus
@w
and so F 0 is nonsingular.
On the other hand, the Jacobian being supplied (or approximated) by the user is
Taking before, and recalling the scaling matrix S defined by (2.6), the corresponding
scaling in x is
I a
Thus we find
Comparing (3.10) and (3.12), we again expect -
J to work well as an approximation to F 0 (x) in a
modified Newton iteration to solve F The corrections to x now take the form
where
It is the class of problems given by (3.8) that we take as the scope of Initialization Problem 1,
for which we have implemented the algorithm described above.
3.4. General index-one systems. Note that (3.8) does not include all fully-implicit index-
one DAEs, because the rank and dependency conditions combined exclude certain index-one sys-
tems. A simple example is the system
This has index one, and it is well-posed for any given value of y 1 (t 0 ). But it does not fit into the
scheme of (3.8), because it contains the derivatives of both variables, and the rank of the 2 \Theta 2
matrix H y 0 is only 1.
In principle, this scheme can be applied to more general index-1 DAE systems by allowing P
to be a more general matrix, not just a permutation. If a constant nonsingular matrix P can be
found that transforms y into
full rank, as in (3.8), then G(t; defined by (3.5). The
vector x and function F (x) can be defined by (3.4) and (3.9) as before, in which y 0;d is defined from
the input vector y 0 as P d y 0 , as in (3.6). The Jacobian F 0 and the approximation -
again
satisfy (3.10) and (3.12). Thus our procedure will still work, as long as the problem supplied by
the user is well-posed. However, in contrast to the case where P is a permutation, once a solution
vector x is found, the vector y may differ from y 0;in in all of its components.
For the example system above, an appropriate matrix P is
making . With this choice, we have
and P
If the input initial value of y 2 differs from g 2 (t 0 ), both components of y 0 would be altered by the
procedure. On the other hand, if preserving y 1 approach would have to be
used.
The procedure with more general constant P determines a consistent set of initial conditions
only in the sense that the initial values of the transformed variables are consistent. Initial values
for the derivatives of the original variables may not be consistent. For example, in the system
2 is initialized correctly, whereas y 0
1 or y 0
2 individually can only be determined
if information about the constraint derivative g 0
2 (t) is available [13].
For the most general index-1 systems, Due to the complications and expense
of finding such a smooth P and continuing with this change of variables in later time steps, we
have chosen not to implement this extension of our algorithm, instead restricting P to the class of
permutation matrices. If necessary, the user may be able to bring the problem to the form (3.8)
by a change of variables as described above.
4. Convergence Theory. In the case of Initialization Problem 1, given by (3.8)-(3.9), the
question of convergence of the Newton or modified Newton iteration arises. Here we give a convergence
analysis for the iteration. The theorem below includes both full and modified Newton
iteration for the same problem, where "full" refers to the fact that the approximate Jacobian is
evaluated at every iteration. We first prove the main convergence result (Theorem 4.1 below),
and then discuss its applicability to the initialization problems of Sections 2 and 3. In the last
subsection, we comment on the use of Newton-Krylov iteration to solve the above initialization
problems.
4.1. Newton iteration convergence. We consider the convergence of the x iteration
and we include both types of Newton iterations by taking
ae
x k for full Newton iteration, or
- x for modified Newton iteration
x is fixed. The function F (x) is a general function here, and it is assumed that there exist
smooth matrix-valued functions -
J(x; h) and C(x) such that -
h of interest.
The following theorem says that this iteration converges under mild smoothness assumptions
on the functions F and C. The norm k \Delta k used here is arbitrary.
Theorem 4.1. Let F : R N ! R N be continuously differentiable in an open convex set D ae R N .
Assume that
(a) there exists x   2 D such that F (x
(c) there exists ae ? 0 such that the neighborhood N(x
condition in N(x   ; ae) with constant fl; and
(d) there exist matrix-valued functions -
J(x; h) and C(x), for all h - 0 and x 2 D, related by
with C(x) satisfying a Lipschitz condition in N(x   ; ae) with constant fl c .
Then there exist constants ffl ? 0 and - h ? 0 such that for x h, and
any - x 2 N(x   ; ffl) in the modified Newton case, the sequence generated by (4.1) is well-defined and
converges to x   . Under these conditions, the iterates obey
)k. In the full Newton case, and in the modified
Newton case,
The proof is an extension of the proof in Dennis and Schnabel [11], p. 90, which treats the
full Newton case with exact Jacobian. The main complication is due to the inaccuracy in the
iteration matrix. Before giving the proof, we give three lemmas that will be useful. In all three,
the hypotheses of the theorem are assumed.
Lemma 4.2. If ffl - minfae; 1=(2fifl)g, -
Proof. First, note that by the Lipschitz condition on F 0 ,
Then, by the perturbation relation (3.1.20) in Dennis and Schnabel [11], p. 45, F 0 (-x) is nonsingular
and
Lemma 4.3. If ffl - minfae; 1=(2fifl)g, then for any -
Proof. By Lemma 4.2, we have using (4.3)
we have
Again, by the perturbation relation (3.1.20) in [11], -
J(-x; h) is nonsingular and
Lemma 4.4. For ffl - minfae; 1=(2fifl)g,
Proof. We have
The norm of the first term is bounded by
xk using Lemma 4.3 and the Lipschitz condition
on F 0 . The second term is bounded by 4hfiC 0 since kC(-x)k - C 0 as in the previous lemma.
Proof of Theorem 4.1. We will derive the recurrence (4.4), which gives linear convergence in
all cases, and q-superlinear convergence in the limit as
Let
ae
oe
Then by Lemmas 4.2 and 4.3, F 0 (x) and -
J(x; h) are both nonsingular with
h. The value of ffl is also a measure of the
nonlinearity of the problem.
We will prove the bound (4.4) by induction on the iteration number k. For 0 - h - h and any
well-defined and we have
Of the two terms in the final equation (4.6), the first is bounded using Lemma 4.1.12 of Dennis and
Schnabel [11], p. 75, which says that kF (x
the norm of the first term is bounded by 2fiflkx  . The second term can be bounded using
Lemma 4.4, which gives a bound of 4fi(flkx
on the norm of the matrix coefficient.
Thus we obtain
and we have shown (4.4) for
Note that since - x 0 and x 0 are both in N(x   ; ffl), we have
using (4.5), we have
and so the Jacobian -
well-defined and nonsingular.
The general induction step goes exactly as with
we obtain
which is (4.4). Again, we have
Thus x showing that the iteration is well-defined. The fact that
gives the linear convergence of the sequence fx k g. 2
4.2. Application to DAE systems. As explained in Section 2, our central objective is to
solve both of the initial condition problems with the help of mechanisms already in place for the
integration of the DAE system itself. The general time step within the DASSL and DASPK solvers
involves the solution of linear systems J \Deltax = R in which the matrix is the DAE system iteration
In order to use values of J inside a Newton or
modified Newton iteration for solving the nonlinear problem F must derive a relationship
between the Jacobian matrix F 0 (x) and J . We have the following lemma, which generalizes similar
relationships derived in earlier sections.
Lemma 4.5. Assume there exists a projection matrix P d such that
G(t;
Define the projection matrix P a j I \Gamma P d and the function F (x) by
with y 0;d given. We then have
for all x.
Proof. From (4.11), we have
Next,
with all the partial derivatives evaluated at
all (t; This then proves the lemma.
In DASPK, starting from input initial guesses y 0 and y 0
corresponding
initial value of x is x
We pick a suitably small value of h, set 1=h, and define
J(x; h) to be the matrix in Lemma 4.5, which satisfies
(4.
with
Note that -
so in the Newton iteration (4.1) based on -
J(x; h), the
correction is
In DASPK, when direct methods are selected, then J(t; supplied by the user
(possibly in approximate form), or generated by difference quotients, and J \Gamma1 is realized by the LU
method. In this case,
0 ) is fixed, and therefore (4.14) represents a modified Newton
method Theorem 4.1 can be applied to this iteration. If one assumes differentiability of
G with respect to y and y 0 , and that the partial derivatives of G with respect to y and y 0 are locally
Lipschitz continuous, then for well-posed initialization problems of the type discussed in Sections 2
and 3, it is clear that the assumptions of the theorem hold. Thus, the iteration on x will converge
given h small enough and a good enough initial guess x 0 .
4.3. Newton-Krylov iteration. When using a Newton-Krylov iteration to solve the initialization
problems of Sections 2 and 3, we use preconditioned GMRES as the linear iteration with
finite-difference approximations involving G(t; to approximate the action of J on an arbitrary
vector, and the preconditioner approximates J \Gamma1 . Once GMRES computes an approximate
solution p k such that small enough, the step \Deltax k is given by \Deltax
This then implies that
small. The reason finite differences of G are used above
(instead of F ) is that the GMRES solver in DASPK uses finite differences of G to approximate J
times an arbitrary vector for the DAE time step. Thus, we are able to apply machinery that is
already available, and the user only need be concerned with J , not F 0 . Because of the relationship
between -
J and F 0 given by Lemma 4.5, it is also clear that a good preconditioner for J will suffice
in the iteration.
The above considerations lead us to consider the convergence of the following inexact Newton
iteration for a general function F (x), where for
with 1. In the general inexact Newton setting, the manner in which the step increment
s k is computed is unimportant. We only need to know that such an s k can be found. Once
again, we assume there exist smooth matrix-valued functions -
J(x; h) and C(x) such that -
of interest.
The following theorem says that this iteration converges under mild smoothness assumptions
on the functions F and C.
Theorem 4.6. Let F : R N ! R N be continuously differentiable in an open convex set D ae R N .
Assume that conditions (a)-(d) of Theorem 4.1 hold. Then there exist constants ffl ? 0 and -
such that for x
h, and any 0 - 1, the sequence generated by (4.15)
is well-defined and converges linearly to x   .
Proof. Since the assumptions (a)-(d) of Theorem 4.1 hold, we can choose ffl and - h as in (4.5).
Thus, for x h, we have -
is nonsingular, the existence of an s 0 satisfying
guaranteed. Given such an s 0 , we have
ks
Hence, it follows that
ks
1, we can choose h smaller (if necessary) so that j +4fi(1+j)hC 0 ! 1. With the residual
norm condition (4.16) on s 0 , the rest of the proof follows along the lines of the proof of Theorem 2.3
in Dembo et. al. [10].
5. The Linesearch Algorithm. In order to improve the robustness of the Newton algorithm
discussed above, we also employ a linesearch backtracking algorithm. Consider a general function
Let x and ffi in R N be such that the residual norm condition
holds, where is the Euclidean norm. The vector ffi can be thought of as an
approximate solution of the Newton equations F 0 obtained using either a direct
solve with an approximate Jacobian or an iterative method such as GMRES. Given such a ffi , it is
shown in Brown and Saad [6] that ffi is a descent direction for f at x, i.e., that
Given ffi a descent direction for f at x, we employ the following backtracking algorithm.
Algorithm 5.1: Given
1.
2. If f(x exit. Otherwise, go to the
next step.
3. Choose - 2 [' min -; ' max -], set -
-, and go to step 2.
The global convergence of this algorithm used in connection with an inexact Newton iteration is
discussed at length in [6]. The simplest choice for the ''s is to take ' 1=2, and this
makes -=2. We use this choice in the implementation here, and also use so that
only a small decrease in f is required.
Given f defined in (5.1), we have and so in the above algorithm
if ffi is the exact Newton step at x, i.e.,
Hence, in this case
and the condition in Step 2 of Algorithm 5.1 is simply
Next, if is the GMRES solution at the m-th step when applied to F 0
then it is shown in [5] that
and the condition in Step 2 of Algorithm 5.1 is
In the present context, we want to solve the nonlinear system F by (2.9) or (3.9).
However, we have no direct measure on the size of the F that is directly available. The weighted
root-mean-square (WRMS) norm used in DASPK for norms of y suggest that we solve instead the
problem
where the matrix A is the current approximate system Jacobian matrix J of (2.4) in the direct
case, or the preconditioner P in the Krylov case, and D is a diagonal matrix containing the weights
to be used in the WRMS norm. Thus, ~
F is just the square of the WRMS norm of A \Gamma1 F . In
either case, it is likely that ~
F (x) is well-scaled in the WRMS norm.
The direction vector ffi that is available to us is the Newton correction given by (2.10) or (3.13).
Thus
S from (3.12) in the case of Initialization Problem 1, or -
in the case of Initialization Problem 2. In both cases, we expect -
but the question arises
then as to whether or not ffi will be a descent direction for ~
F at the current approximate solution
x. An easy calculation gives
using ~
if we can assure that -
I is small, it follows that ffi will
be a descent direction for ~
F at x.
6. Implementation. We implemented the algorithms described above for Initialization Problems
1 and 2 as new options in the general-purpose DAE solver DASPK [7]. Initialization Problem
1 has been implemented for the more general class of index-one systems described in Section 3.3.
Here we will give a few details concerning the implementation, and describe briefly how to use the
new options.
We will assume here that the reader is familiar with the use of DASSL and DASPK. Detailed
descriptions of those solvers can be found in [2] and [7], respectively. These solvers normally require
the initial t, y, y 0 to be consistent. Within the dependent variable vector Y in DASPK, we denote
by Y d the differential variables, and denote by Y a the algebraic variables. By specifying the input
parameter INFO(11), DASPK will solve one of the following two initialization problems:
Initial values are already consistent (default).
Given Y d , calculate Y a and Y 0
d . If this option
is specified, the user must identify for DASPK the differential and algebraic components
of Y. This is done by setting (for I
is a differential variable, and
is an algebraic variable.
2. Solve Initialization Problem 2: Given Y 0 , calculate Y .
In either case, initial values for the given components are input, and initial guesses for the unknown
components must also be provided as input.
The algorithm for Problem 1 requires an initial stepsize or scaling, h, to determine
As a first approximation, we try the initial stepsize h 0 which is used by DASSL and DASPK ([2],
p. 128). Since we do not know in advance if this value is small enough to achieve convergence of
the modified Newton iteration, we set up a loop, starting with . If the initialization fails,
we divide h by 10. If the initialization fails for MXNH (nominally = 5) different values of h, the
code returns an error flag to the user program. In our experience, if the initialization succeeds, it
usually succeeds with the initial choice h 0 . For Initialization Problem 2, we always set
there is no such loop on h. In either case, once the initialization has been completed, we reset the
initial stepsize h 0 for the first step of DASPK based on the newly computed initial values, using
the formula in [2].
For a given value of c, the initialization problem is solved with either a modified Newton
method or an inexact Newton method [10] similar to that used in the general time step. It has
been augmented by the linesearch algorithm described above for improved global convergence. In
both the case of Initialization Problem 1 as given in Eq. (3.9), and Problem 2 as given by Eq.
(2.9), we must solve a system
where F (x) is the residual of the DAE system at t 0 , y, y 0 , and x represents the variables for which we
are solving. A Jacobian matrix (or preconditioner matrix, in the Krylov case) is obtained either by
finite-difference approximations or by calling a user subroutine to calculate the Jacobian, depending
on the option specified. The routines to specify the DAE and the Jacobian/preconditioner are
exactly the same ones which are needed for the time integration. The Newton iteration is given in
terms of an approximation to the system Jacobian J by
where - is the relaxation steplength (0 ! - 1) from the linesearch algorithm given in Section
5. (For Initialization Problem 2, the scaling matrix -
S is absent.) The code is organized so that
the Newton solver is independent of which initialization problem is being solved. After a vector
has been calculated, a separate routine is called to construct and apply the increment
\Deltax as follows:
ffl For Problem 1, we increment y 0 by \Gamma-P a p and increment y 0
0 by \Gamma-cP d p.
ffl For Problem 2, we increment y 0 by
The meaning of J in (6.2) depends on the choice of methods. In the case of direct methods,
the value of J is the approximation to the system Jacobian evaluated (by difference quotients or
user-supplied routine) at the start of the iteration. Then the evaluation of J involves a call
to a back-substitution routine. In the case of Krylov methods, however, J refers to the exact system
Jacobian at the current values of y and y 0 , and J \Gamma1 F (x) is evaluated by a call to a routine for the
preconditioned GMRES method [16]. This makes use of the preconditioner P supplied by the user,
evaluated at the start of the iteration. In both cases, the inaccuracy of the fixed approximation (J
or P ) to the system Jacobian is an additional potential cause of difficulty for the Newton iteration.
The complete algorithm actually involves three loop levels for Problem 1, and two levels for
Problem 2. At the innermost level, up to MXNIT Newton iterations are performed with a given
value of h and a given value of the Jacobian or preconditioner. The iteration is considered to have
converged if the scaled residual is small in norm:
ae
where A is the current approximate system Jacobian J in the direct case, and the preconditioner
matrix P in the Krylov case. Here the test constant EPCONI is
is the tolerance for the Newton iteration in the subsequent time steps, and
EPINIT is a "swing factor" nominally equal to 0.01. The norm used throughout is the weighted
root-mean-square (WRMS) norm in which the weights are formed from the user-supplied tolerances
(see [7]).
The values ae m from the mth iteration are used to infer a convergence rate,
If convergence is not achieved in MXNIT iterations, the strategy for repeated
attempts depends on RATE. If convergence failed, but RATE - 0:8 (the iterations are converging,
but slowly), then we retry the Newton iteration with the current values of y and y 0 and a new value
for A (i.e. for J or P ), up to a limit of MXNJ such attempts. In addition, in the case of the Krylov
method, if the GMRES solver failed to converge after at least two Newton iterations, but RATE
1, the Newton iteration is retried with a new value for P . If the limit of MXNJ retries is reached,
we reduce h and retry the iteration (again with a new A and the current y and y 0 ) in the case of
Problem 1, or give up and return an error flag in the case of Problem 2. If convergence of the inner
Newton iteration failed but RATE ? 0:8 (or some other recoverable failure occurred), we retry the
iteration with a reduced value of h and the initial y and y 0 (Problem 1) or give up (Problem 2).
The total number of iterations performed can therefore be as large as MXNH*MXNJ*MXNIT in
Problem 1, and MXNJ*MXNIT in Problem 2.
Currently we have set We have set
the case of direct methods, and in the case of Krylov methods. However,
all four of these controls are optional inputs to DASPK, so that a user may specify different values.
In addition, an option is provided to turn off the linesearch algorithm.
Actually, an additional level of logic has been added around the initialization algorithm described
above. The reason for it is that the error weights involved in all convergence and error tests
depend on the current solution vector:
Thus, while the initialization algorithm may have succeeded using the weights evaluated at the
initial guess, those weights may differ greatly from the updated values using the converged y
vector. Since updating the weights at every iteration seems rather extreme, we have adopted the
following scheme: With weights set using the input y vector, the initialization algorithm is called,
and if it succeeds, we update the weights and call it a second time. If it again succeeds, we update
the weights again, and proceed to the first time step. If either initialization fails, an error flag
is returned to the user. In the case of the Krylov method, on the second initialization call, the
preconditioner is not updated unless and until there is a convergence failure.
7. Numerical Experiments. We tested the initialization algorithm on several problems and
found that it performed much as expected. In the course of development and debugging, we used a
simple index-one system of size 2, having a known analytic solution. For both the first and second
initialization problem types, and for a wide range of initial guesses, the initialization algorithm
converged within the limits imposed, for both the direct and Krylov method options. All attempts
to integrate the system without the initialization option failed except when the initial values were
consistent.
For a more realistic test, we used a model of a multi-species food web [3], in which mutual
competition and/or predator-prey relationships in a spatial domain are simulated. Here we consider
a 2-species model, species 1 being the prey and species 2 being the predator, and with the
predator assumed to have an infinitely fast reaction rate. Specifically, the model equations for the
concentration vector are:
with
The interaction and diffusion coefficients (a could be functions of (x; general. The
choices made for this test problem are as follows:
and
The domain is the unit square 10. The boundary conditions are of
homogeneous Neumann type (zero normal derivatives) everywhere. The coefficients are such that
a unique stable equilibrium is guaranteed to exist when derivatives appear in
the equations for species 2 [3]. Empirically, a stable equilibrium appears to exist for (7.1) when ff
and fi are positive, although it may not be unique. In our tests on this problem we take
and 100, for which there is considerable spatial variation in the solution.
The PDE system (7.1), together with the boundary conditions, was discretized with central
differencing on an L \Theta L mesh, as described in [7]. We have taken which is quite sufficient
for accurate spatial resolution. The resulting DAE system G(t; Y; Y 0
The tolerances used were were run on a Sun Sparc-10 workstation.
7.1. Initialization Problem 1. In the tests on this problem reported in [7], the initial conditions
were taken to be mildly peaked functions that nearly satisfy the constraint equations:
ae
The predator value c determined by the equation f 2 (x; is an approximate quasi-steady
state (QSS) value. The original DASPK solver has no difficulty with this problem, without
further adjustment of the initial values. However, we expect that in a typical application of this
type it is impractical to find such accurate initial values. So for our tests, we will prescribe a flat
value
predas the initial guess in the input Y array, and invoke the new algorithm for Initialization Problem
1. For the present problem parameters, the QSS values of c 2 at time are all within 10% of
so we vary c pred
We will report here only tests with the Krylov method (GMRES) option in DASPK, and as
a preconditioner we use a product of a spatially-based factor and a reaction-based factor. In the
notation of [7], this is given by
Here R and S are (respectively) the reaction and diffusion terms of the right-hand side of the
DAE system, so that the problem has the form G(t; Y; Y 0
I
I 1 is the identity matrix
with 0 in place of 1 in positions corresponding to the components c 2 ). The spatial factor in P SR
consists of 5 Gauss-Seidel iterations, and the reaction factor uses difference quotient approximations
for the diagonal blocks. For the DASPK input parameters relating to the Krylov method, default
values were specified.
In

Table

7.1 below, we summarize the results of the DASPK tests with the new initialization
algorithm incorporated in it. For each value of c pred
(with QSS denoting the values in (7.6)), the
tabulated quantities are:
number of Newton iterations in the initial condition calculation
number of linear iterations in the initial condition calculation
total number of Newton iterations to complete the integration
total number of linear iterations to complete the integration
total number of residual evaluations to complete the integration.
The numbers NNI0 and NLI0 measure the cost of the initialization algorithm, while NNI, NLI,
and NRE measure the total cost of solving the problem. Convergence (to correct values) in the
initialization was achieved at a very reasonable additional cost for pred
Evidently,
the convergence region for the initialization of this problem is strongly skewed to the high side, but
does permit errors of at least 40% on the low side. In the case c pred
the algorithm converged,
but to the value c which corresponds to a solution that is valid but different from the one of
interest here.
For comparison, when the initial condition calculation option was not selected, only the QSS
initial values were successful, and in that case the total cost figures were
709. These are slightly larger than with the initialization, indicating that even the
approximate QSS values from (7.6) are somewhat in error. The unmodified DASPK solver, when
run with its initial condition option on, was also unable to solve any case except the QSS initial
values, and in that case the total costs were 971. The failed cases
either halted in the initialization algorithm, or (when the initialization option was off) failed in the
first time step with either repeated corrector convergence failures or repeated error test failures.
2 On the basis of experience with these tests, however, we have changed the default value of NRMAX, the maximum
number of GMRES restarts, from 2 to 5.
c pred
failed in I.C. calculation

Table
Test results for new initialization algorithm on food web problem
7.2. Initialization Problem 2. In these tests, we specify the initial time derivatives y 0 to be
0, i.e., we are posing the steady-state problem for (7.1). Since we have no explicit time-dependence
in the right-hand sides, once the consistent initial values are determined, the solution to the DAE
problem is constant in time. This type of problem is much harder than Initialization Problem 1,
because of the absence of the time derivative operator. Moreover, the preconditioners devised in
for the DAE problem itself are less useful here.
We first describe tests using the direct method. These specified a banded Jacobian, generated
internally by difference quotients, where the two half-bandwidths are equal to
simplicity, the initial guesses for the discrete c i values were taken to be spatially flat values with
prey
Because the subsequent time integration is not an issue here, we stopped it at We
performed tests for a variety of values of ff and fi, revealing, as in the case of Problem 1, a nontrivial
region of convergence in each case. Table 7.2 below (upper half) gives the results for the case
where the tabulated counter NNI0 is defined as before. Convergence is achieved
with no difficulty (always using the full Newton step) for (at least) the values prey
In
an interval about c prey
and in the interval c prey
the algorithm fails to find a solution. In
an interval about c prey
converges to an incorrect solution (that has negative values of c i ).
In an interval about c prey
converges to the correct solution, but with difficulty, in that the
linesearch algorithm must choose vales of - ! 1. For reference, we note that the true steady state
values of c 1 in this case range from 9.9 to 66.
For this 2-D problem, using a Jacobian with the full bandwidth is quite costly. In an attempt
to reduce costs, we also tested with half-bandwidths equal to 1, corresponding to an approximate
Jacobian that ignores the diffusion terms. However, the results were completely unsuccessful. The
resulting lumped tridiagonal preconditioner is evidently too inaccurate.
In considering tests with the Krylov method for this problem, the choice of a preconditioner
is problematical. In terms of the form
I the true Jacobian for the steady state
problem is Since the initialization algorithm sets user-supplied
preconditioner, the choice P SR of (7.7) used for Problem 1 is undefined. We therefore use
\GammaR Y , a block-diagonal matrix involving only the reaction Jacobian elements. We again tried a
variety of values of ff and fi, and provided flat initial guesses (7.8). However, for the larger values
of these parameters, it was found that convergence of the GMRES iteration was much slower than
in the case of Initialization Problem 1. This is to be expected, since the diffusion terms contribute
significantly to the system but are completely absent in the preconditioner. In order to achieve
convergence, we therefore increased the Krylov subspace parameters over their default values,
setting the maximum size of the Krylov subspace (MAXL) to 20, and the number of GMRES
restarts allowed (NRMAX) to 19. This allows a total of 400 GMRES iterations on each linear
system. For values of c prey
the algorithm appears to fail, while for values of prey
(at least), it converges to the correct solution.
c prey
linear method NNI0 NLI0 Notes
direct - fails in I.C. calculation
direct 11 - incorrect solution
direct 11 - linesearch min
direct - fails in I.C. calculation
50 direct 11 -
direct 11 -
direct 12 -
50 Krylov - fails in I.C. calculation

Table
Test results for food web problem, Initialization Problem 2
In all of the cases tabulated, we compared the computed solution vector from the initialization
algorithm with that from a more accurate solution with the direct method, integrated with tighter
tolerances to is virtually at steady state. All of the values from the Problem 2
tests had errors less than the tolerances imposed. For example, for prey
the maximum relative error observed was about
Users of DASPK should be cautioned that Initialization Problem 2 is more difficult than Initialization
Problem 1, and that some extra effort may be necessary. By comparison with Problem
1 and with the time integration, convergence of the algorithm for Problem 2 is more sensitive to
the initial guess and to the quality of the approximate Jacobian J or preconditioner P . If the J
or P used in the time steps is a good approximation only in the limit c ! 1 as was the
case in the Problem 1 food web tests, a different preconditioner for the steady-state initialization
problem (where should be seriously considered. The user can easily determine in JAC and
PSOL whether the preconditioner has been called for a steady-state initial condition calculation
because the parameter CJ will be equal to zero (it is nonzero in any other situation) and branch
accordingly to the appropriate preconditioner.



--R

Developing Software for Time Dependent Problems using the Method of Lines and Differential-Algebraic Integrators
Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations
Decay to Uniform States in Food Webs
A Local Convergence Theory for Combined Inexact-Newton/Finite-Difference Projection Methods
Hybrid Krylov Methods for Nonlinear Systems of Equations
Convergence Theory of Nonlinear Newton-Krylov Algorithms
Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems
Consistent Initial Conditions for Linear Time Varying Singular Systems
A Computational Method for General Higher-Index Nonlinear Singular Systems of Differential Equations

Numerical Methods for Unconstrained Optimization and Nonlinear Equations
A Shooting Method for Fully-Implicit Index-2 Differential-Algebraic Equations
Approximation Methods for the Consistent Initialization of Differential-Algebraic Equations
Zur Theory und Numerischen Realisierung von L-osungmethoden bei Differentialgleichungen mit Angekoppelten Algebraischen Gleichung
The Consistent Initialization of Differential-Algebraic Systems
GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems
--TR

--CTR
Alan C. Hindmarsh , Peter N. Brown , Keith E. Grant , Steven L. Lee , Radu Serban , Dan E. Shumaker , Carol S. Woodward, SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers, ACM Transactions on Mathematical Software (TOMS), v.31 n.3, p.363-396, September 2005
C. T. H. Baker , C. A. H. Paul , H. Tian, Differential algebraic equations with after-effect, Journal of Computational and Applied Mathematics, v.140 n.1-2, p.63-80, 1 March 2002
D. A. van Beek , V. Bos , J. E. Rooda, Declaration of unknowns in DAE-based hybrid system specification, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.13 n.1, p.39-61, January
differential-algebraic equations;consistent initial conditions