Mesh Independence of Matrix-Free Methods for Path Following.
In this paper we consider a matrix-free path following algorithm for nonlinear parameter-dependent compact fixed point problems. We show that if these problems are discretized so that certain collective compactness and strong convergence properties hold, then this algorithm can follow smooth folds and capture simple bifurcations in a mesh-independent way.
Introduction
. The purpose of this paper is to extend the results in [4], [5], and [16] on
mesh-independent convergence of the GMRES, [21], iterative method for linear equations to a
class of matrix-free methods for solution of parameter dependent nonlinear equations of the form
In Banach space.
We present an algorithm for numerical path following and detection of simple bifurcations
together with conditions on a sequence of approximate problems,
(with G consistency of notation) that imply that the performance of the algorithm is independent
of the level h of the discretization. Hence, for such problems, methods, and discretizations,
the difficulties raised in [23] and [24] will not arise.
In this section we set the notation and specify the kinds of singlarities that we will consider.
This discussion, and that of algorithms for path following in x 2.1 and detection of simple bifurcation
and branch switching in x 3, do not depend on the discretization and we use the notation G
for both G 0 and G h for h ? 0. When we describe our assumptions on the discretization in x 4 and
present an example in x 5, the distinction between G 0 and G h becomes important and the notation
in those sections reflects that difference.
1.1. Notation and Simple Singularities. We let
denote the solution path in X \ThetaIR of pairs (u; -) that satisfy (1.1), G u denotes the Fr-echet derivative
with respect to u, and G - the derivative with respect to the scalar -.
Version of May 19, 1998.
y Department of Applied Mathematics, National Chiao-Tung University, Hsin-Chu 30050, Taiwan.
(ferng@math.nctu.edu.tw). The research of this author was supported by National Science Council, Taiwan.
z 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-9700569.
Throughout this paper we make
ASSUMPTION 1.1. G is continuously Fr- echet differentiable in (u; -). G u is a Fredholm
operator of index zero, and
Here COM denotes the space of compact operators.
For A 2 L(X), the space of bounded operators on X , we let N (A) be the null space of A and
R(A) the range of A. Assumption 1.1 implies that R(G u (u; -)) is closed in X and the dimension
of N (G u (u; -)) is the co-dimension of R(G u (u; -)).
Following [14], we say that a point (u;
ffl is a regular point if G u (u; -) is nonsingular,
ffl is a simple singularity if 0 is an eigenvalue of G u with algebraic and geometric multiplicity
one,
- is a simple fold if it is a simple singularity and G - 62 R(G u ).
- is a simple bifurcation point if it is a simple singularity and G - 2 R(G u ).
2. Algorithms for Path Following and Branch Switching. In this section we describe the
iterative methods for path following, detection of bifurcation, and branch switching that we analyze
in the subsequent sections and discuss some alternative approaches. Our approach is matrix-free,
which means that we use no matrix storage or matrix factorizations at all. However, there are
many related methods to which our mesh-independence analysis is applicable. We refer the reader
to [13] and [7] for other methods for path following and branch switching and to [20], [25], and
[10] for iterative methods for detection of bifurcation.
2.1. Path Following and Arclength Continuation. For path following in the absence of singularities
the standard approach is a predictor-corrector method. The methods differ in the manner
in which u 0 Assume that solution point and we want
to compute at a nearby - 1 .
If the Jacobian G u nonsingular, then the implicit function theorem insures the existence
of a unique smooth arc of solutions (u(-), through
more, with the smoothness assumption on G, it follows that u 0
d-
u(-) exists and differentiation
of (1.1) with respect to - gives the equation for u 0 .
G u
The Euler predictor uses u 0 directly and approximates
While the results in this paper can be applied to the solution of (2.1) by GMRES, we choose
to avoid the addtional solve and approximate u 0 with a difference. The secant predictor is
solution pair with -
CONTINUATION METHODS 3
The nonlinear iteration is an inexact Newton iteration [6] in which
where d k satisfies the inexact Newton condition
kG
The inexact Newton condition can be viewed as a relative residual termination criterion for an
iterative linear solver applied to the equation for the Newton step
G
When GMRES is used as the linear iteration, as it is in this paper, the nonlinear solver is referred
to as Newton-GMRES.
The same procedure can be used if there are simple folds if the problem is expanded by using
pseudo-arclength continuation, [12], [13]. Here we introduce a new parameter s and solve
is a normalization equation.
We will use the secant normalization
when two points on the path are available and the norm-based normalization
to begin the path following. Both normalizations are independent of the discretization. This independence
plays a role in the analyis that follows.
It is known, [14], that arclength continuation turns simple folds in (u; -) to regular points in
3. Singular Point Detection and Branch Switching. Although a continuation procedure incorporated
with a pseudo-arclength normalization can circumvent the computational difficulties
caused by turning points and usually jump over bifurcation points, it is usually desirable and important
to be able to detect and locate the singularities in both cases. In the case where direct
methods are used for linear algebra a sign change in the determinant of F can be used to detect
simple bifurcation and a change in sign of d-=ds to detect simple folds. We use a variation
of the approach in [10], which requires only the largest (in magnitude) solution of a generalized
eigenvalue problem.
Suppose that (x(s a ); s a ) and (x(s b ); s b ) are two regular points and the path following procedure
is going from s a to s b . Let A(s a
A(s a ) and A(s b ) are nonsingular. Recall that if (x(s); s) is a bifurcation point on the solution path,
then the null space N dimensional and there exists a nonzero vector v such that
Applying the Lagrange interpolation to
A(s a
where the matrix E(s) denotes the perturbation matrix in the interpolation. Then instead of solving
the usually nonlinear eigenvalue problem (3.1) we make a linear approximation
A(s a
which leads to a generalized eigenvalue problem
A(s a
with
Therefore, an s value which causes the Jacobian F x (x(s); s) to be singular can be approximated
by
Our approach differs from that in [10] in that the roles of s a and s b are interchanged. The
method of [10] solves
where
In the approach of [10], F x (x; s) is factored at s a in order to compute dx=ds for the predictor; that
factorization is used again to precondition an Arnoldi iteration for the corrector, and one last time
in the formulation of the eigenvalue problem to detect singularities. Since we do not factor F x at
all, the roles of s a and s b can be interchanged. We found in our experiments that our approach,
using (3.3), gave a better approximation of the location of the singularity that one using (3.6).
Equation (3.4) implies that a bifurcation point s closest to the interval [s a ; s b ] corresponds to
the largest eigenvalue in magnitude oe of (3.3). A negative oe signals that there is a bifurcation
point between s a and s b , indicates that there is a bifurcation point close behind s a ,
and a "large" positive oe means that a bifurcation point is approaching. When oe - 1, it should be
interpreted that no bifurcation point is nearby.
regular point and A(s b ) is nonsingular, the generalized eigenvalue problem
(3.3) can be solved via the equivalent linear eigenvalue problem
(3.
CONTINUATION METHODS 5
where only the largest eigenvalue and corresponding eigenvector are needed. We will solve (3.8)
with the Arnoldi method and prove that if s a and s b are not too near a singlar point then the
eigenvalue problem is well conditioned in a mesh independent way.
A simple fold point can be predicted in a similar way with s replaced by - and A(s) replaced
by
At a simple bifurcation point two branches of solutions intersect nontangentially. We next
describe how the information obtained from the solution of the eigenvalue problem can be used for
branch switching. Suppose a bifurcation point x on the primary branch
is determined. Since the eigenvector w corresponding to the largest eigenvalue oe for the eigenvalue
problem described above is an approximation for a null vector of F x
On the other hand, the tangent vector -
which is the solution of
G
is also a null vector of [G u
x 0 are
linearly independent it is recommended in [10], where G u is factored and dx=ds is computed
using (3.9), that N ([G u approximated by spanf -
wg. In the matrix-free case
considered in this paper, we approximate the tangent vector by a secant approximation
We approximate the tangent direction of the new branch by orthogonalizing w against ffi x
To obtain a regular point on the secondary branch, we solve the following augmented nonlinear
system
where
and ffl is a "switching factor". This switching factor is problem dependent and should be chosen
large enough so that the solution of (3.11) does not fall back to the primary branch [14]. The
nonlinear equation (3.11) can be solved by Newton-type method with initial iterate x
After
moving onto the secondary branch, the continuation procedure for path following on the secondary
branch is identical to the one on the primary branch.
In [10] the linear systems for computing tangent vectors (3.9) are solved with a preconditioned
Arnoldi iteration and the eigenvalues of the Hessenberg matrix produced in the Arnoldi iteration
are used to predict singular points. Since a factorization of A(s a ) is used as the preconditioner
in [10], the prediction comes with very little cost. In our approach, Jacobians are not factored
and the tangent vector is not computed, instead, a secant vector is used as an approximation.
The eigenvalue problem for singularity prediction is performed separately using some iterative
algorithm. This leaves us the flexibility of choosing any robust iterative solver and the associated
preconditioner for the continuation procedure.
Other methods for singularity detection based on solution of eigenvalue problems have been
proposed in [25], [9], and [20].
4. Approximations and Mesh Independence. This is the only section in which the properties
of the discrete problems are explicitly addressed. We assume, as is standard in the integral
equations literature [1], [3], that the discretization used in (1.2) has been constructed, by interpolation
if necessary, so that G h has the same domain and range as G.
Recall [1] that a family of operators fT ff g ff2A is collectively compact if the set
is precompact in X . In (4.1) B is the unit ball in X . In the special case that the indexing set A is an
interval [0; h 0 ], we will denote the index by h. We say that strongly as h ! 0 and write
in the norm of X for all u 2 X .
All of the results are based on the following results from [1] and some simple consequences.
THEOREM 4.1. Let fT h g be a family of collectively compact operators on X that converge
strongly to T 2 COM(X). Assume that I \Gamma T is nonsingular. Then there is h 0 ? 0 such that
nonsingular for all h - h 0 and converges strongly to
A simple compactness argument implies uniform bounds for parameter dependent families of
collectively compact strongly convergent operators.
COROLLARY 4.2. Let a ! b be given and assume that
is a collectively compact family of operators such that T h
for each fixed s.
Assume that fT h (s)g is uniformly Lipschitz continuous in s and that I \Gamma T (s) is nonsingular for
all s 2 [a; b]. Then there are M; h 0 ? 0 such that
for all h - h 0 and s 2 [a; b].
We let
We assume that
ASSUMPTION 4.1.
1. G h (u; -) ! G(u; -) for all (u; -) 2 X \Theta IR as h ! 0.
2. G h is Lipschitz continuously Fr- echet differentiable and the Lipschitz constants of G h
u and
are independent of h.
3. For all (u; -) 2 X \Theta IR
(4.
CONTINUATION METHODS 7
4. There are
then the families of operators
are collectively compact.
We augment G with a mesh-independent Lipschitz continuously differentiable arclength
normalization defining
0:
Examples of normalizations N that do not depend on h are (2.6) and (2.7). We let
Consistently with the notation in the previous sections, we let u h (-) and x h solutions
to G h (u;
4.1. The Corrector Equation and Simple Folds. At regular points when - is used as the
continuation parameter the equation for the Newton step is
The theory in [1] asserts that if G u (u; -) is nonsingular, so is G h
sufficiently small h and
strongly convergent to G u (u; -) \Gamma1 . However, this convergence is not uniform and in
order to invoke the results of [1] and [5] we must take care to remain away from singular points.
One can treat simple folds as regular points by means of arclength continuation. If we solve the
augmented system, (2.5), and this system has only regular points, then
x
k is bounded on finite
segments of \Gamma. If, moreover, the normalization N(u; -; s) does not depend on the discretization, as
it will not if (2.6) or (2.7) are used, then our assumptions imply that the finite dimensional problems
are as well conditioned as the infinite dimensional problems.
THEOREM 4.3. Let Assumptions 1.1 and 4.1 hold. Let [s a ; s b ] be such that
is a single smooth arc and F x is nonsingular on - \Gamma(s) then there are ffi 1 , K, and h 0 such that
Proof. Assumption 1.1 and the mesh independence of N imply that we may apply Corollary
4.2 to T h
x (x(s); s). Hence, there is such that for all h - h 0 and
x
Let L be the (h-independent) Lipschitz constant of F h
x
. Then if
Moreover
and hence if M Lffi ! 1 the Banach Lemma implies that F h
is nonsingular and
1=(2ML) the proof is complete with
The bound, (4.7), part 2 of Assumption 4.1, and the Kantorovich theorem [11], [15], [17],
imply convergence of x h to x.
COROLLARY 4.4. Let Assumptions 1.1 and 4.1 hold. Let [s a ; s b ] be such that f -
[s a ; s b ]g is a single smooth arc with at most simple fold singularities. Then there is a unique
solution arc for F h ,
and x h
So, for h sufficiently small, the secant predictor
converges to x (0) uniformly in (s Hence if the steps in arclength fffi n
s
are independent of h and the secant predictor is used, then the accuracy of the initial iterate to the
corrector equation and by Theorem 4.3 the condition of the linear equation for the Newton step is
independent of h. Hence the performance of the nonlinear iteration is independent of h.
As for the GMRES iteration that computes the Newton step, the methods from [16], [4],
and [5], may be extended to show that the GMRES iteration for the Newton step converges r-
superlinearly in a manner that is independent of h, x, and s. All that one needs is an uniform
clustering of the eigenvalues of F h
on M(ffi) that is independent of x, s, and h. This follows
directly from a theorem in [1]. In the results that follow, we count eigenvalues by multiplicity and
order them by decreasing distance from 0.
THEOREM 4.5. Let fT h g be a family of collectively compact operators on X that converge
strongly to K 2 COM(X). Let f- h
be the eigenvalues of T h and f- j g 1
be the eigenvalues
of T . Let ae ? 0 and assume that . Then there is
then
. Moreover, if - 1 has algebraic and geometric multiplicity one, then so does
. Moreover there is a sequence of eigenfunctions fw h g of T h corresponding to the
eigenvalue - h
1 so that w h ! w, and eigenfunction of T correpsonding to the eigenvalue - 1 .
With this in hand, the eigenvalue clustering result can be obtained in the same way that Theorem
4.3 was derived from Theorem 4.1. We begin with the analog of Corollary 4.2.
COROLLARY 4.6. Let the assumptions of Theorems 4.3 and 4.5 hold. Let fmu h
be
the eigenvalues, counted by multiplicity, of I \Gamma F x (x; s). Let ae ? 0 be given. There are h
such that if h - h 0 then
CONTINUATION METHODS 9
THEOREM 4.7. Let the assumptions of Theorems 4.3 and 4.5 hold. Let fmu h
be the
eigenvalues, counted by multiplicity, of I \Gamma F x (x; s). Let ae ? 0 be given. There are
such that if h - h 0 and (x; s) 2 N
for all
Before we state our r-superlinear convergence theorem we must set some more notation. We
write the equation for the Newton step for the corrector equation as
We let d h
denote the kth GMRES iteration and r h
the kth GMRES residual. As is
standard in nonlinear equations, d h
s). The proof of Theorem 4.8 is
similar to that of the similar results in [4].
THEOREM 4.8. Let Assumptions 1.1 and 4.1 hold. Let [s a ; s b ] be such that -
[s a ; s b ]g is a single smooth arc and F x is nonsingular on -
\Gamma(s) Then there are
a continuous function M on (0; 1) such that for all ae 2 (0; 1) and (x; s) ae N (ffi 3 ) \Theta [s a ; s b ] and
Theorem 4.8 states that any desired rate of linear convergence can be obtained in a mesh
independent way and hence, for (x; s) fixed, the convergence is r-superlinear in a mesh independent
way.
4.2. Simple Bifurcation. Now consider a path -
which has a single simple bifurcation at be the branch of solutions
that intersects -
\Gamma 0 at (x(s c ); s c ). We will denote solutions on - \Gamma 1 by y and the arclength parameter
on -
t. Hence, for
We set (y(t c ); t c ) be the bifurcation point on Hence
For any ffi ? 0, the results in x 4.1 hold for the paths
We define
0 in a similar way. Since x h
the results in the previous section, simple bifurcation from - \Gamma hcan only arise for s 2
In this section we describe how the prediction of the bifurcation point, the conditioning of the
generalized eigenproblem, the solution of that eigenproblem, and the accurate tracking of the other
branch depend on h. Since perturbation of G, even in finite dimension, can change the structure
of -
\Gamma 0 from two intersecting arcs to two disconnected arcs [18], [19], we cannot show that the
bifurcation will be preserved, however if s c is far enough away from s a and s b , the new other path
will be detected even if it does not correspond to a bifurcation for the finite dimensional problem.
For
will be nonsingular and we can consider the eigenvalue problems
where the largest eigenvalue in magnitude is sought. If (4.11) is solved via the Arnoldi method, the
matrix-vector products of F h
with a vector will each require a product
of F h
e. a linear solve. If this solve is performed with GMRES,
then Theorem 4.8 implies that the number of iterations required to approximate that matrix-vector
product can be bounded independently of h; x; s.
Our assumptions imply, if s b \Gamma s a is sufficiently small, that oe = oe 0 is a simple eigenvalue with
geometric multiplicity one. Since
is a collectively compact and strongly convergent sequence, we can apply Theorem 4.5 to conclude
that for h sufficiently small, oe h is a simple eigenvalue with geometric multiplicity one as well.
Moreover oe h will be well separated from the next largest eigenvalue in magnitude, hence the
conditioning of the eigenvalue problem is independent of h.
5. Numerical Results. In the example presented in this section we use the normalization
equation
are two previously computed solutions and
used in the experiments. We used a step of
in s for path following.
We apply the secant predictor (2.3) to generate the initial iterate. For the corrector, we choose
a forward-difference Newton-GMRES algorithm [15]. The outer iteration that generates the sequence
when the nonlinear residual norm satisfies
where the absolute error tolerance - were used. To
avoid oversolving on the linear equation for the Newton step z k the forcing terms j k such that
CONTINUATION METHODS 11
were adjusted with a method from [8]. We use the l 2 norm and scale the norm by a factor of
1=N for differential and integral equations so that the results are independent of the computational
mesh.
The generalized eigenvalue problem (3.3) that characterizes the singularity prediction and
branch switching need not be solved at each new point on the path is computed (i. e. s
is too frequent). In the example considered in this section, we solved the eigenvalue problem to
make a predicition of bifurcation after a step of
had been taken. We used a simple version of the Arnoldi, [2], method based on modified Gram-Schmidt
orthogonalization to solve (3.3) with reorthogonalization at each iteration. The Jacobian
of the inflated system is usually neither symmetric nor positive definite, we treat the generalized
eigenvalue problem as linear eigenvalue problem and solve the linear system involved explicitly
with preconditioned GMRES.
After j steps the Arnoldi method produces
is an upper Hessenberg matrix with the h ij 's as its nonzero entries and Q j is orthogonal
with the q j 's as columns. If ('; y) is an eigenpair of H j and y, then
provides a computable error bound. The Ritz pair, ('; w), is used to approximate the eigenpair of
We terminated the Arnoldi iteration when
5. The prediction -
- is accepted when - lies in between - a and - b , or
We found this strategy sufficient for the example here and were able to use (3.10) to move onto the
new branch. Of course, the values of \Delta eig is problem dependent (as it is for other methods) as is
the criteria for accepting a prediction.
As an example we consider the following two-point boundary value problem [14].
Uniform mesh of spacing
is used to define the grid points ft j g N
as
order finite difference discretization leads to the discrete nonlinear system of order
FIG. 5.1. Solution u( 1) versus -.
Solution
at the boundaries.
Starting from the trivial solution the algorithm traces the solution path, locates
bifurcation points, switches branches, and captures the secondary solution curves. Indeed there are
total two bifurcation points at - 81; \Gamma81 and eight turning points at - \Sigma11; \Sigma110; \Sigma336
in the region of interest. The primary solution branch represents the symmetric solutions and
the secondary solution branch represents nonsymmetric periodic solutions bifurcating from the
primary branch. Figure 5.1 plots the solution u( showing the folds and bifurcations on
the curve.
In

Table

5.1 we report the bifurcation predictions along the solution path going toward -
\Gamma81. One can see that (3.3) is a more accurate predictor than (3.6). Tables 5.2 and 5.3 illustrate
the mesh independence of the linear and nonlinear iterations at various points on the primary (5.2)
and secondary ((5.3) branches. The iteration statistics remain virtually unchanged as the mesh is
refined.

Table

5.4 lists the total number of preconditioned GMRES iterations and Arnoldi iterations
corresponding to three representative prediction intervals. In Table 5.4 we show how both the
GMRES iterations needed to approximate the product of F h
u with G u in the case where we predict a turning point) and the overall number of Arnoldi
iterations is independent of the mesh. Each Arnoldi step requires about 12-14 GMRES iterations.
This is similar to the numbers listed in Table 5.2 and 5.3. The residual in the table is the Arnoldi
residual when the iteration terminates.
CONTINUATION METHODS 13

Prediction of bifurcation point along the primary branch going toward - \Gamma81.
- a - b prediction (3.3) prediction (3.6)
6.0998e+00 1.0574e+01 1.2614e+02 -9.3188e+03
8.4208e+00 1.1361e+00 -2.2855e+01 -3.5346e+01
1.1361e+00 -7.1699e+00 -3.6599e+01 -4.7426e+01

Total number of Newton and preconditioned GMRES iterations on the primary branch.
problem size - value Newton P-GMRES

Total number of Newton and preconditioned GMRES iterations at switching point and on the secondary branch.
problem size - value Newton P-GMRES
switching 7 21
switching 7 22
switching 6 20
Total number of preconditioned GMRES and Arnoldi iterations required at different sections on the path when a
prediction procedure is performed.
path N - a - b P-GMRES Arnoldi residual
turning 64 7.0814e+00 1.0632e+01 53 5 2.9490e-02
point 128 6.0998e+00 1.0574e+01 50 5 1.9156e-05
near-by 256 9.1379e+00 1.0883e+00
regular 64 -3.1304e+01 -4.0215e+01
point 128 -2.4500e+01 -3.3295e+01 69 5 3.3738e-05
bifurcation
point 128 -6.8734e+01 -7.7621e+01 56 4 1.3263e-06
detected 256 -7.4823e+01 -8.3605e+01 62 4 1.3251e-06

Acknowledgments

. The authors wish to thank Gene Allgower for making us aware of reference
[22], Klaus B-ohmer for providing us with a copy, and Alastair Spence for many other pointers
to the literature. We are also grateful to Dan Sorensen for his council on the solution of eigenvalue
problems.



--R

Collectively Compact Operator Approximation Theory
The principle of minimized iterations in the solution the matrix eigenvalue problem
A survey of numerical methods for Fredholm integral equations of the second kind
GMRES and the minimal polynomial
Convergence estimates for solution of integral equations with GMRES

Lecture Notes on Numerical Analysis of Bifurcation Problems.
Choosing the forcing terms in an inexact Newton method
WANG, Nonequivalence deflation for the solution of matrix latent value problems
A new algorithm for numerical path following applied to an example from hydrodynamic flow
Functional Analysis
Constructive methods for bifurcation and nonlinear eigenvalue problems

Lectures on Numerical Methods in Bifurcation Theory
Iterative methods for linear and nonlinear equations
GMRES and integral operators
Iterative Solution of Nonlinear Equations in Several Variables
Column buckling-an elementary example of bifurcation

Algorithms for the nonlinear eigenvalue problem
GMRES a generalized minimal residual algorithm for solving nonsymmetric linear systems
Andwendung von Krylov-Verfahren auf Verzweigungs-und Fortsetzungsprobleme
An adaption of Krylov subspace methods to path following problems

Numerical path following and eigenvalue criteria for branch switching
--TR
fold point;bifurcation;arnoldi method;singularity;matrix-free method;GMRES;path following;collective compactness;mesh independence