Convergence Analysis of Pseudo-Transient Continuation.
Pseudo-transient continuation ($\Psi$tc) is a well-known and physically motivated technique for computation of steady state solutions of time-dependent partial differential equations. Standard globalization strategies such as line search or trust region methods often stagnate at local minima. \ptc succeeds in many of these cases by taking advantage of the underlying PDE structure of the problem. Though widely employed, the convergence of \ptc is rarely discussed. In this paper we prove convergence for a generic form of \ptc and illustrate it with two practical strategies.
Introduction
. Pseudo-transient continuation (\Psitc) is a method for computation of steady-state
solutions of partial differential equations. We shall interpret the method in the context of
a method-of-lines solution, in which the equation is discretized in space and the resulting finite
dimensional system of ordinary differential equations is written as x
@x
@t
and the the discretized spatial derivatives are contained in the nonlinear term F (x). Marching out
the steady state by following the physical transient may be unnecessarily time consuming if the
intermediate states are not of interest. On the other hand, Newton's method for F
usually not suffice, as initial iterates sufficiently near the root are usually not available. Standard
globalization strategies [12, 17, 25], such as line search or trust region methods often stagnate at
local minima of kFk [20]. This is particularly the case when the solution has complex features
such as shocks that are not present in the initial iterate (see [24], for example). \Psitc succeeds in
many of these cases by taking advantage of the PDE structure of the problem.
1.1. The Basic Algorithm. In the simple form considered in this paper, \Psitc numerically integrates
the initial value problem
to steady state using a variable timestep scheme that attempts to increase the timestep as F (x)
approaches 0. V is a nonsingular matrix used to improve the scaling of the problem. It is typically
diagonal and approximately equilibrates the local CFL number (based on local cell diameter and
local wave speed) throughout the domain. In a multicomponent system of PDEs, not already
Version of April 20, 1997.
y North Carolina State University, Department of Mathematics and Center for Research in Scientific Computation,
Box 8205, Raleigh, N. C. 27695-8205 (Tim Kelley@ncsu.edu). The research of this author was supported by
National Science Foundation grant #DMS-9321938.
z Computer Science Department, Old Dominion University, Norfolk, VA 23519-0162 (keyes@cs.odu.edu)
and ICASE, MS 132C, NASA Langley Research Center, Hampton, Virginia 23681-0001. The research of this author
was supported by National Science Foundation grant #ECS-9527169, and NASA grants NAGI-1692 and NAS1-19480,
the latter while the author was in residence at ICASE.
properly nondimensionalized, V might be a block diagonal matrix, with blocksize equal to the
number of components.
We can define the \Psitc sequence fx n g by
where F 0 is the Jacobian (or the Fr-echet derivative in the infinite-dimensional situation). Algorithmically

ALGORITHM 1.1.
1.
2. While kF (x)k is too large.
(a) Solve (ffi
(b)
(c) Evaluate F (x).
(d) Update ffi.
The linear equation
for the Newton step that is solved in step 2a is the discretization of an PDE and is usually very large.
As such, it is typically solved by an iterative method which terminates on small linear residuals.
This results in an inexact method [11], and a convergence theory for \Psitc must account for this; we
provide such a theory in x 3. We have not yet explained how ffi is updated, nor explored the reuse
of Jacobian information from previous iterations. These options must be considered to explain the
performance of \Psitc as practiced, and we take them up in x 3, too. In order to most clearly explain
the ideas, however, we begin our analysis in x 2 with the most simple variant of \Psitc possible,
solving (1.3) exactly, and then extend those results to the more general situation.
As a method for integration in time, \Psitc is a Rosenbrock method ([14], p. 223) if ffi is fixed.
One may also think of this as a predictor-corrector method, where the simple predictor (result from
the previous timestep) and a Newton corrector are used. To see this consider the implicit Euler step
from x n with timestep
z
In this formulation z n+1 would be the root of
Finding a root of G with Newton's method would lead to the iteration
If we take - the first Newton iterate is
This leaves open the possibility of taking more corrector iterations, which would lead to a different
form of \Psitc than that given by (1.2). This may improve stability for some problems [16].
PSEUDO-TRANSIENT CONTINUATION 3
1.2. Time Step Control. We assume that ffi n is computed by some formula like the "switched
evolution relaxation" (SER) method, so named in [21], and used in, e.g. [19], [24], and [33]. In
its simplest, unprotected form, SER increases the timestep in inverse proportion to the residual
reduction.
Relation (1.5) implies that, for n - 1,
In some work [16], ffi n is kept below a large, finite bound ffi max . Sometimes ffi n is set to 1
(called "switchover to steady-state form" in [13]) when the computed value of
In view of these practices, we will allow for somewhat more generality in the formulation of the
sequence fffi n g. We will assume that ffi 0 is given and that
for n - 1. The choice in [24] and [33] (equation (1.5)) is
OE(-:
Other choices could either limit the growth of ffi or allow ffi to become infinite after finitely many
steps. Our formal assumption on OE accounts for all of these possibilities.
ASSUMPTION 1.1.
1.
2. Either -
and the timesteps are held bounded. If - then switchover to steady-state form
is permitted after a finite number of timesteps.
In [16] the timesteps are based not on the norms of the nonlinear residuals kF (x n )k but on
the norms of the steps This policy has risks in that small steps need not imply small
residuals or accurate results. However if the Jacobians are uniformly well conditioned, then small
steps do imply that the iteration is near a root. Here formulae of the type
are used, where OE satisfies Assumption 1.1.
1.3. Iteration Phases. We divide the \Psitc iteration into three conceptually different and separately
addressed phases.
1. The initial phase. Here ffi is small and x is far from a steady state solution. This phase
is analyzed in x 2.3. Success in this phase is governed by stability and accuracy of the
temporal integration scheme and proper choice of initial data.
2. The midrange. This begins with an accurate solution x and a timestep ffi that may be small
and produces an accurate x and a large ffi. We analyze this in x 2.2. To allow ffi to grow without
loss of accuracy in x we make a linear stability assumption (part 3 of Assumption 2.1).
3. The terminal phase. Here ffi is large and x is near a steady state solution. This is a small
part of the overall process, usually requiring only a few iterations. Aside from the attention
that must be paid to the updating rules for ffi, the iteration is a variation of the chord method
[25, 17].
We analyze the terminal phase first, as it is the easiest, in x 2.1. Unlike the other two phases,
the analysis of the terminal phase does not depend on the dynamics of x (x). The
initial and midrange phases are considered in x 2.3, with the midrange phase considered first to
motivate the goals of the initial phase. This decomposition is similar to that proposed for GMRES
and related iterations in [22] and is supported by the residual norm plots reported in [24, 10].
2. Exact Newton Iteration. In this section we analyze the three phases of the solver in reverse
order. This ordering makes it clear how the output of an earlier phase must conform to the demands
of the later phase.
2.1. Local Convergence: Terminal Phase. The terminal phase of the iteration can be analyzed
without use of the differential equation at all.
LEMMA 2.1. Let fffi n g be given by either (1.6) or (1.8) and let Assumption 1.1 hold. Let
be nonsingular, and F 0 be Lipschitz continuous with Lipschitz constant fl in a
ball of radius ffl about x   .
Then there are then the
sequence defined by (1.2) and (1.6) satisfies
and x
Proof. Let denote the error. As is standard, [12], [17], we analyze convergence in
terms of the transition from a current iterate x c to a new iterate x+ . We must also keep track of the
change in ffi and will let ffi c and ffi + be the current and new pseudo-timesteps.
The standard analysis of the chord method ([17], p. 76) implies that there are ffl 1 - ffl and K C
such that if
The constant K C depends only ffl 1 , F , and x   and does not increase if ffl 1 is reduced.
Now let \Delta \Gamma1and ffl 1 be small enough to satisfy
and, in particular, F increase
if needed so that
PSEUDO-TRANSIENT CONTINUATION 5
where - denotes condition number. Equations (2.1) and (2.2) imply that
If fffi n g is computed with (1.6) we use the following inequality from [17] (p. 72)
and (2.3) to obtain
2:
We then have by Assumption 1.1 that
t is from Assumption 1.1.
If fffi n g is computed with (1.8), we note that
and hence
as before.
In either case, Therefore we may continue the iteration
and conclude that at least q-linearly with q-factor of 1=2.
If we complete the proof by observing that since
superlinear convergence follows from (2.1).
The following simple corollary of (2.1) applies to the choice OE(-.
COROLLARY 2.2. Let the assumptions of Lemma 2.1 hold. Assume that OE(-.
1. Then the convergence of fx n g to x   is q-quadratic.
2.2. Increasing the Time Step: Midrange Phase. Lemma 2.1 states that the second phase
of \Psitc should produce both a large ffi and an accurate solution. We show how this happens if the
initial phase can provide only an accurate solution. This clarifies the role of the second phase
in increasing the timestep. We do this by showing that if the steady state problem has a stable
solution, then \Psitc behaves well. We now make assumptions that not only involve the nonlinear
function but also the initial data and the dynamics of the IVP (1.1).
ASSUMPTION 2.1.
1. F is everywhere defined and Lipschitz continuously Fr- echet differentiable, and there is
2. There is a root x   of F at which F 0
then the solution of the initial value problem
converges to x   as t !1.
3. There are ffl 2
for all
The analysis of the midrange uses part 3 of Assumption 2.1 in an important way to guarantee
stability. The method for updating ffi is not important for this result.
THEOREM 2.3. Let fffi n g be given by either (1.6) or (1.8) and let Assumption 1.1 hold. Let
Assumption 2.1 hold. Let ffi max be large enough for the conclusions of Lemma 2.1 to hold. Then
there is an ffl 3 ? 0 such that if
or x
Proof. Let 1 is from Lemma 2.1 and ffl 2 is from part 3 of Assumption
2.1. Note that
Now there is a c ? 0 such that
for all x such that kek - ffl 1 . Hence, reducing ffl 3 further if needed so that ffl 3 ! fi=(2c), we have
If
for all n - 1 and hence x n converges to x   q-linearly with q-factor (1
This convergence implies that x
(1.8) is used.
This result says that once the iteration has found a sufficiently good solution, either convergence
will happen or the the iteration will stagnate with inf latter failure mode is, of
course, easy to detect. Moreover, the radius ffl 3 of the ball about the root in Theorem 2.3 does not
depend on inf ffi n .
2.3. Integration to Steady State: Initial Phase. Theorem 2.3 requires an accurate estimate
of x   , but asks nothing of the timestep. In this section we show that if ffi 0 is sufficiently small, and
(1.6) is used to update the timestep, then the dynamics of (1.1) will be tracked sufficiently well for
such an approximate solution to be found. It is not clear how (1.8) allows for this.
THEOREM 2.4. Let fffi n g be given by (1.6) and let Assumption 1.1 hold. Let Assumption 2.1
hold. There is a -
ffi such that if ffi 0 - ffi then there is an n such that
Proof. Let S be the trajectory of the solution to (2.5). By Assumption 2.1 x   satisfies the
assumptions [12, 17, 25], for local quadratic convergence of Newton's method and therefore there
are ffl 4 and ffl f such that if
then suffice for the conclusions of Theorem 2.3 to hold. Let
We will show that if ffi 0 sufficiently small, then . By
(1.7), if
then
as long as kF (x)k - ffl f and x k is within ffl 4 of the trajectory S.
Let z be the solution of (2.5). Let be such that for all t ? T ,
Consider the approximate integration of (2.5) by (1.2). Set
If holds. This cannot happen if n ?
(which implies that t n ? T ). Therefore the proof will be complete if we can show that
for all
Note that (1.2) may be written as
There is an m 1 such that the last term in (2.9) satisfies, for sufficiently small and
)k. Then we have, by our assumptions on F , that there is an
such that
Finally, there is an m 2 such that for for sufficiently small and
Setting we have for all n - 1 (as long as ffi n is sufficiently small and
As long as (2.7) holds, this implies that
Consequently, as is standard [14, 15],
and using (2.8),
then This completes the proof.
The problem with application of this proof to the update given by (1.8) is that bounds on ffi like
(2.7) do not follow from the update formula.
3. Inexact Newton Iteration. In this section we look at \Psitc as implemented in practice. There
are two significant differences between the simple version in x 2 and realistic implementations:
1. The Fr-echet derivative
recomputed with every timestep.
2. The equation for the Newton step is solved only inexactly.
Before showing how the results in x 2 are affected by these differences, we provide some more
detail and motivation.
Item 1 is subtle. If one is solving the equation for the Newton step with a direct method,
then evaluation and factorization of the Jacobian matrix is not done at every timestep. This is a
common feature of many ODE and DAE codes, [30, 26, 27, 3]. Jacobian updating is an issue in
continuation methods [31, 28], and implementations of the chord and Shamanskii [29] methods for
general nonlinear equations [2, 17, 25]. When the Jacobian is slowly varying as a function of time
or the continuation parameter, sporadic updating of the Jacobian leads to significant performance
gains. One must decide when to evaluate and factor the Jacobian using iteration statistics and (in
the ODE and DAE case) estimates of truncation error. Temporal truncation error is not of interest
to us, of course, if we seek only the steady-state solution.
CONTINUATION 9
In [16] a Jacobian corresponding to a lower-order discretization than that for the residual was
used in the early phases of the iteration and in [19], in the context of a matrix-free Newton method,
the same was used as a preconditioner.
The risks in the use of inaccurate Jacobian information are that termination decisions for
the Newton iteration and the decision to reevaluate and refactor the Jacobian are related and one
can be misled by rapidly varying and ill-conditioned Jacobians into premature termination of the
nonlinear iteration [30, 32, 18]. In the case of iterative methods, item 1 should be interpreted to
mean that preconditioning information (such as an incomplete factorization) is not computed at
every timestep.
means that the equation for the Newton step is solved inexactly in the sense of [11], so
that instead of
where s is given by (1.3), step s satisfies
for some small j, which may change as the iteration progresses. Item 1 can also be viewed as
an inexact Newton method with j reflecting the difference between the approximate and exact
Jacobians.
The theory in x 2 is not changed much if inexact computation of the step is allowed. The proof
of Lemma 2.1 is affected in (2.1), which must be changed to
This changes the statement of the lemma to
LEMMA 3.1. Let fffi n g be given by either (1.6) or (1.8) and let Assumption 1.1. Let F (x
be nonsingular, and F 0 be Lipschitz continuous with Lipschitz constant fl in a ball of radius
ffl about x   .
Then there are ffl 1 ? 0, -
j for all n, and
then the sequence defined by (3.1), (3.2), and (1.6) satisfies
and x
Corollary 2.2 becomes
COROLLARY 3.2. Let the assumptions of Lemma 3.1 hold. Assume that OE(-.
1. Then the convergence of fx n g to x   is q-superlinear if
locally q-quadratic if
The analysis of the midrange phase changes in (2.6), where we obtain
for some K j ? 0. This means that -
j must be small enough to maintain the q-linear convergence
of fx n g during this phase. The inexact form of Theorem 2.3 is
THEOREM 3.3. Let fx n g be given by (3.1) and (3.2) and let fffi n g be given by either (1.6)
or (1.8). Let Assumption 1.1 hold. Let Assumption 2.1 hold. Let ffi max be large enough for the
conclusions of Lemma 2.1 to hold. Then there are ffl 3 ? 0 and - j such that if j n -
and
or x
Inexact Newton methods, in particular Newton-Krylov solvers, have been applied to ODE/DAE
solvers in [1], [5], [4], [6], [7], and [9]. The context here is different in that the nonlinear residual
F (x) does not reflect the error in the transient solution but in the steady state solution.
The analysis of the initial phase changes through (2.10). We must now estimate
and hence, assuming that the operators are uniformly bounded, there is m 3
such that
ks
and hence
We express (3.5) as
Hence, if and the inexact form of Theorem
2.4:
THEOREM 3.4. Let fx n g be given by (3.1) and (3.2) and let fffi n g be given by (1.6) and let
Assumption 1.1 hold. Let Assumption 2.1 hold. Assume that the operators
are uniformly bounded in n. Let ffl ? 0. There are - ffi and -
j such that if
j then there
is an n such that
The restrictions on j in Theorem 3.4 seem to be stronger than those on the results on the
midrange and terminal phases. This is consistent with the tight defaults on the forcing terms for
methods when applied in the ODE/DAE context [1, 5, 6, 7, 9].
4. Numerical Experiments. In this section we examine a commonly used \Psitc technique,
switched evolution/relaxation (SER) [21], applied to a Newton-like method for inviscid compressible
flow over a four-element airfoil in two dimensions. Three phases corresponding roughly to
the theoretically-motivated iteration phases of x 2 may be identified. We also compare SER with a
different \Psitc technique based on bounding temporal truncation error (TTE) [20]. TTE is slightly
PSEUDO-TRANSIENT CONTINUATION 11
-0.4
Zoomed Grid
FIG. 4.1. Unstructured grid around four-element airfoil in landing configuration - near-field view.
more aggressive than SER in building up the time step in this particular problem, but the behavior
of the system is qualitatively the same.
The physical problem, its discretization, and its algorithmic treatment in both a nonlinear
defect correction iteration and in a Newton-Krylov-Schwarz iteration - as well as its suitability
for parallel computation - have been documented in earlier papers, most recently [10] and the
references therein. Our description is correspondingly compact.
The unknowns of the problem are nodal values of the fluid density, velocities, and specific
total energy, at N vertices in an unstructured grid of triangular cells
(see Fig. 4.1). The system F discretization of the steady Euler equations:
r
r
r
where the pressure p is supplied from the ideal gas law, and fl is the
ratio of specific heats. The discretization is based on a control volume approach, in which the
control volumes are the duals of the triangular cells - nonoverlapping polygons surrounding each
vertex whose perimeter segments join adjacent cell centers to midpoints of incident cell edges.
Integrals of (4.1)-(4.3) over the control volumes are transformed by the divergence theorem to
contour integrals of fluxes, which are estimated numerically through an upwind scheme of Roe
type. The effective scaling matrix V for the \Psitc term is a diagonal matrix that depends upon the
mesh.
The boundary conditions correspond to landing configuration conditions: subsonic (Mach
number of 0.2) with a high angle of attack of (5 ffi ). The full adaptively clustered unstructured grid
contains 6,019 vertices, with four degrees of freedom per vertex (giving 24,076 as the algebraic
dimension of the discrete nonlinear problem). Figure 4.1 shows only a near-field zoom on the
full grid, whose far-field boundaries are approximately twenty chords away. The initial pseudo-
\Gamma4 corresponds to a CFL number of 20. The pseudo-timestep is allowed to
grow up to six orders of magnitude over the course of the iterations. It is ultimately bounded at
guaranteeing a modest diagonal contribution that aids the invertibility of (ffi
The initial iterate is a uniform flow, based on the far field boundary conditions - constant
density and energy, and constant velocity at a given angle of attack.
The solution algorithm is a hybrid between a nonlinear defect correction and a full Newton
method, a distinction which requires further discussion of the processes that supply F (x) and F 0 (x)
within the code. The form of the vector-valued function F (x) determines the quality of the solution
and is always discretized to required accuracy (second-order in this paper). The form of the
approximate Jacobian matrix F 0 (x), together with the scaling matrix V and time step ffi, determines
the rate at which the solution is achieved but does not affect the quality of a converged result, and
is therefore a candidate for replacement with a matrix that is more convenient. In practice, we perform
the matrix inversion in (1.2) by Krylov iteration, which requires only the action of F 0 (x) on
a series of Krylov vectors, and not the actual elements of F 0 (x). The Krylov method was restarted
preconditioned with 1-cell overlap additive Schwarz (8 subdomains).
Following [5, 8], we use matrix-free Fr-echet approximations of the required action:
However, when preconditioning the solution of (1.2), we use a more economical matrix than the
Jacobian based on the true F (x), obtained from a first-order discretization of the governing Euler
system. This not only decreases the number of elements in the preconditioner, relative to a true
Jacobian, but also the computation and (in the parallel context) communication in applying the
preconditioner. It also results in a more numerically diffusive and stable matrix, which is desirable
for inversion. The price for these advantages is that the preconditioning is inconsistent with the
true Jacobian, so more linear subiterations may be required to meet a given linear convergence
tolerance. This has an indirect feedback on the nonlinear convergence rate, since we limit the work
performed in any linear subiteration in an inexact Newton sense.
In previous work on the Euler and Navier-Stokes equations [10, 23], we have noted that a \Psitc
method based on a consistent high-order Jacobian stumbles around with a nonmonotonic steady-state
residual norm at the outset of the nonlinear iterations for a typical convenient initial iterate far
from the steady-state solution. On the other hand, a simple defect correction approach, in which
is based on a regularizing first-order discretization everywhere it appears in the solution of
(1.2), not just in the preconditioning, allows the residual to drop smoothly from the outset. In
this work, we employ a hybrid strategy, in which defect correction is used until the residual norm
has fallen by three orders of magnitude, and inexact Newton thereafter. As noted in x 3, inexact
FIG. 4.2. SER convergence history
iteration based on the true Jacobian and iteration with an inconsistent Jacobian can both be gathered
under the j of (3.2), so the theory extends in principal to both.
With this background we turn our attention to Fig. 4.2, in which are plotted on a logarithmic
scale against the \Psitc iteration number: the steady-state residual norm jjF (x n )jj 2 at the beginning
of each iteration, the norm of the update vector jjx and the pseudo-timestep ffi n .
The residual norm falls nearly monotonically, as does the norm of the solution update. Asymptotic
convergence cannot be expected to be quadratic or superlinear, since we do not enforce
in (3.5). However, linear convergence is steep, and our experience shows that overall
execution time is increased if too many linear iterations are employed in order to enforce
asymptotically. In the results shown in this section, the inner linear convergence tolerance was set
at 10 \Gamma2 for the defect correction part of the trajectory, and at 10 \Gamma3 for the Newton part. The work
was also limited to a maximum of 12 restart cycles of 20 Krylov vectors each.
Examination of the pseudo-timestep history shows monotonic growth that is gradual through
the defect correction phase (ending at rapidly growing, and asymptotically
at (beginning at show momentary retreats from ffi max in
response to a refinement on the \Psitc strategy that automatically cuts back the pseudo-timestep by
a fixed factor if a nonlinear residual reduction of less than 3is achieved at the exhaustion of the
maximum number of Krylov restarts in the previous step (during the terminal Newton phase).
Close examination usually reveals a stagnation plateau in the linear solver, and it is more cost
effective to fall back to the physical transient to proceed than to grind on the ill-conditioned linear
problem. These glitches in the convergence of jjF are not of nonlinear origin.
Another timestep policy, common in the ODE literature, is based on controlling temporal
truncation error estimates. Though we do not need to maintain temporal truncation errors at low
levels when we are not attempting to follow physical transients, we may maintain them at high
levels as a heuristic form of stepsize control. This policy seems rare in external aerodynamic
simulations, but is popular in the combustion community and is implemented in [20]. The first
FIG. 4.3. TTE convergence history
neglected term in the Euler discretization of @x
@t
is reasonable mixed absolute-relative
bound on the error in the i th component of x at the n th step is
where
can be approximated
Taking - as 3and implementing this strategy in the Euler code in place of SER yields the results in
Fig. 4.3. Arrival at ffi max occurs at the same step as for SER, and arrival at the threshold jjF
occurs one iteration earlier. However, the convergence difficulties after having arrived at
are slightly greater.
5. Conclusions. Though the numerical experiments of the previous section do not confirm
the theory in detail, in the sense that we do not verify the estimates in the hypotheses, a reassuring
similarity exists between the observations of the numerics and the conceptual framework of the
theory, which was originally motivated by similar observations in the literature. There is a fairly
long induction phase, in which the initial iterate is guided towards the Newton convergence domain
by remaining close to the physical transient, with relatively small timesteps. There is a terminal
phase which can be made as rapid as the capability of the linear solver permits (which varies from
application to application), in which an iterate in the Newton convergence domain is polished.
Connecting the two is a phase of moderate length during which the time step is built up towards
the Newton limit of ffi max , starting from a reasonably accurate iterate. The division between these
phases is not always clear cut, though exogenous experience suggests that it becomes more so
when the corrector of x 1 is iterated towards convergence on each time step. We plan to examine
PSEUDO-TRANSIENT CONTINUATION 15
this region of parameter space in conjunction with an extension of the theory to mixed steady/\Psitc
systems (analogs of differential-algebraic systems in the ODE context) in the future.

Acknowledgments

. This paper began with a conversation at the DOE Workshop on Iterative
Methods for Large Scale Nonlinear Problems held in Logan, Utah, in the Fall of 1995. The authors
appreciate the efforts and stimulus of the organizers, Homer Walker and Michael Pernice. They
also wish to thank Peter Brown, Rob Nance, and Dinesh Kaushik for several helpful discussions on
this paper. This paper was significantly improved by the comments of a thoughtful and thorough
referee.



--R

The Numerical Solution of Initial Value Problems in Differential-Algebraic Equations
Some efficient algorithms for solving systems of nonlinear equations
VODE: A variable coefficient ODE solver


Using Krylov methods in the solution of large-scale differential-algebraic systems

Hybrid Krylov methods for nonlinear systems of equations
Pragmatic experiments with Krylov methods in the stiff ODE setting


Numerical Methods for Nonlinear Equations and Unconstrained Opti- mization
Towards polyalgorithmic linear system solvers for nonlinear elliptic problems
Numerical Initial Value Problems in Ordinary Differential Equations
Analysis of numerical methods
Robust linear and nonlinear strategies for solution of the transonic Euler equations
Iterative Methods for Linear and Nonlinear Equations
Termination of Newton/chord iterations and the method of lines
Aerodynamic applications of Newton-Krylov-Schwarz solvers
A parallelized elliptic solver for reacting flows
Experiments with implicit upwind methods for the Euler equations
Convergence of Iterations for Linear Equations
Application of Newton-Krylov methodology to a three-dimensional unstructured Euler code
A Newton's method solver for the Navier-Stokes equations
Iterative Solution of Nonlinear Equations in Several Variables
A description of DASSL: a differential/algebraic system solver
Description and use of LSODE
Driven cavity flows by efficient numerical techniques
A modification of Newton's method
Implementation of implicit formulas for the solution of ODEs
An error estimate for the modified newton method with applications to the solution of nonlinear two-point boundary value problems
Accurate and economical solution of the pressure head form of Richards' equation by the method of lines
Newton solution of inviscid and viscous problems
--TR

--CTR
W. K. Anderson , W. D. Gropp , D. K. Kaushik , D. E. Keyes , B. F. Smith, Achieving high sustained performance in an unstructured mesh CFD application, Proceedings of the 1999 ACM/IEEE conference on Supercomputing (CDROM), p.69-es, November 14-19, 1999, Portland, Oregon, United States
T. Coffey , R. J. McMullan , C. T. Kelley , D. S. McRae, Globally convergent algorithms for nonsmooth nonlinear equations in computational fluid dynamics, Journal of Computational and Applied Mathematics, v.152 n.1-2, p.69-81, 1 March
D. Gropp , Dinesh K. Kaushik , David E. Keyes , Barry Smith, Performance modeling and tuning of an unstructured mesh CFD application, Proceedings of the 2000 ACM/IEEE conference on Supercomputing (CDROM), p.34-es, November 04-10, 2000, Dallas, Texas, United States
Feng-Nan Hwang , Xiao-Chuan Cai, A parallel nonlinear additive Schwarz preconditioned inexact Newton algorithm for incompressible Navier-Stokes equations, Journal of Computational Physics, v.204
Howard C. Elman , Victoria E. Howle , John N. Shadid , Ray S. Tuminaro, A parallel block multi-level preconditioner for the 3D incompressible Navier--Stokes equations, Journal of Computational Physics, v.187 n.2, p.504-523, 20 May
Keyes , Lois Curfman Mcinnes , M. D. Tidriri, Globalized Newton-Krylov-Schwarz Algorithms and Software for Parallel Implicit CFD, International Journal of High Performance Computing Applications, v.14 n.2, p.102-136, May       2000
D. A. Knoll , D. E. Keyes, Jacobian-free Newton-Krylov methods: a survey of approaches and applications, Journal of Computational Physics, v.193 n.2, p.357-397, 20 January 2004
pseudo-transient continuation;nonlinear equations;global convergence;steady state solutions