Michael E Henderson  Michael E Henderson photo         

contact information

Principal Research Staff Member, Applied Mathematics
Thomas J. Watson Research Center, Yorktown Heights, NY USA


Professional Associations

Professional Associations:  ACM  |  IEEE   |  Society for Industrial and Applied Mathematics  |  Tau Beta Pi - The Engineering Honor Society

more information

More information:  Multifario  |  Fields Institute lecture (2001) on Multiparameter Continuation

Continuation Methods

Continuation Methods are used to compute solution manifolds of nonlinear systems. They are usually cast in terms of the solution of an equation:

F=0, F:R^n->R^{n-k}

The solutions of these equations are (generically) a curve when k=1, a surface when k=2 and a k-manifold in general. They arise frequently in engineering and scientific problems because those problems are often formulated as determining a function that satisfies some set of equations, for example, the Navier-Stokes equations, Maxwell's equations, or Newton's law.


Argonne has a page describing homotopy methods, with some pointers to available software.

 There are two broad approaches to solving these problems. Both require a starting solution and compute solutions that can be found by continuous deformation. They are simplicial continuation methods and predictor-corrector methods.

 This is a work in progress. It has been fun trying to set all of this down, but it takes a lot of time, and I got more ambitious as I went along.

Continuation Methods

I've got a way of looking at these methods, that I think makes it clearer how the various algorithms work, and why. The idea is that at each step in the continuation a neighborhood A_m(u_m) of a point u_m on the manifold is merged with the previously computed part of the manifold M_m, which produces M_{m+1}. A new point u_{m+1} is then selected from M_{m+1}, and the process is repeated until the part of M inside some finite region is covered. The continuation starts with a single point u_0, provided by the user, and as long as A_m(u_m) has some part which is exterior to M_mthe process will eventually terminate. The natural choice is therefore that u_{m+1} be on or near the boundary of M_m.

A step of a general continuation

The different methods make different choices for how the neighborhoods are computed, and how the point and manifold are represented.

Simplicial Continuation Methods

The system defining M is k equations short of being a square system. Roughly speaking, if k constraints are added, the system has a single point as a solution. So, for example, we might require u to lie on an n-k dimensional linear subspace, or on an n-k dimensional face of a simplex in R^n. This is what simplicial continuation does. It represents u_m as the intersection of M with an n-k dimensional simplex sigma_m. The neighborhood A_m(u_m) is  then the intersection of M and a set of n dimensional simplices with sigma_m in its interior. Of course, the intersection need not be done explicitly, it is enough just to have A_m(u_m) be the n dimensional simplices. If  M_m is represented in the same way, then if the simplices in M_m and A_m(u_m) are compatible, (the intersection of two simplcies is either empty or another simplex), then the merge operation is just adding the simplices from A_m(u_m) which are not already present. To ensure this is the case, M_m and A_m(u_m) are subsets of a reference simplicial decomposition of R^n.

Simplicial continuation in the general context

The next ''point'', sigma_{m+1}, is any n-k dimensional face on the boundary of M_{m+1} which crosses M.
Since an n dimensional simplex has (n+1) choose (n-k+1), n-k dimensional faces, checking all of the n-k dimensional faces of A_m(u_m) for crossing M may require quite a bit of work if n is large (only those exterior to M_m have to be checked).

The algorithm is formulated in terms of a labeling of the vertices of a simplicial complex , and a test of which faces of each simplex are completely labeled. This language obscures a very simple idea. The mapping F(u) assigns a value in R^{n-k} to each vertex of a (n-k)-dimensional simplex in R^n. So for example, a scalar to each endpoint of an edge in the plane-

The image of a triangle containing the origin to a line segment.

It is easy to test if the interval which is created by mapping the edge contains zero, you just test the sign of F at the endpoints. If one is positive and one is negative the edge contains zero.

 In higher dimensions it is the same, but different. If we have an (n-k)-dimensional simplex, S, and F maps from S to R^{n-k}, the question is whether or not the image of S contains the origin. For n=3, and k=1 this is a mapping of one triangle to another --

The image of the faceof  a tetrahedron to a plane

To test if a simplex lies near the origin, there is an integer labeling which assigns to the vertex the number of leading coordinates which are positive (in 2-d points in quadrant (-,+) and (-,-) are labeled "0", points in (+,-) "1" and in (+,+) "2"). If none of the vertices have the same label, |F| must be small within the image. If you create the smallest hypercube containing the image (with faces perpendicular the coordinate axes), then this test is simply that that hypercube has a face on either side of each coordinate axis. With a bound on the size of the original simplex, and a bound of the derivative of F, this yields a bound on the size of the hypercube.

 There is an alternative, called vector labeling. The function values themselves are used as labels, and a linear mapping is defined from a standard simplex,

t_i ge 0,  sum t_i le 1,

to the image of the simplex. This assigns a coordinate to each point inside the image, and if the Jacobian F_u is full rank over the simplex, the mapping can be inverted to find the point which maps to 0.

f =  F(u_0) + sum (F(u_i)-F(u_0))t_i

f(t)=0 is an (n-k)x(n-k) linear system of equations. If the solution t has positive coordinates, and their sum is less than 1, the origin is contained in the image simplex. This approach should be familiar to those of you doing Finite Elements. In 1-d it reduces to finding t in [0,1] so that (1-t) F(u_0)+tF(u_1)=0. Note that we also get an approximation to the zero, from

u = u_0 + sum (u_i-u_0)t_i

There one last issue, and that is that the image of the entire simplex under F has curved edges. If that curvature is too large the origin may be in the simplex formed by the images of the vertices, but not in the image of the original simplex. This requires bounds on the derivatives of F relative to the size of the edges of the simplex.

 Assuming that there is a procedure for testing whether the image of a simplex contains the origin, a continuation can be constructed for a mapping from R^n to R^{n-k} by starting with a given n dimensional simplex, and marking all of the n-k faces which cross the manifold. For n=3, and k=2, which is a surface in three space, we would mark the n-k-cells, or edges, which cross the surface. If we take the set of all simplices in R^n that have been found to contain points on either side of the manifold, and list the n-k faces on the boundary, whose image under F contains zero, we have an approximation to the boundary of the computed part of the manifold. The next step in the continuation takes one of these marked boundary faces, attaches a new simplex to it, and marks the n-k faces that lie on the boundary.


We have so far assumed that the simplicial complex on which the continuation operates is given. Most of us can take a regular grid in 2-d and split each square into two triangles to yield a triangulation of part of the plane. Many of us can take a cubic grid in 3-d and split each cube into five tetrahedra and reflect them across cube faces to give a tetrahedralization of part of 3-space. Generating simplicial grids in more general settings in 2 and 3 dimensions is called Finite Element Grid Generation, and has absorbed many careers. I'd rather not do it in higher dimensions.

Fortunately, there are some standard simplicial complexes available that cover all of R^n (see [4]). The basic operation we need is to find the simplex adjacent to a given simplex across one of it's faces, so it is natural to use triangulations built by reflections across a face (Coxeter enumerated these). However, there are disadvantages to using a fixed decomposition. Real problems have varying scales, and one part of a manifold may be very flat, while others have large curvature. A fixed decomposition requires that the small cells be used in all parts of the manifold (there are ways, but...). Therefore, we might wish to generate a decomposition on the fly. This can be very simple, by choosing a new vertex on a line through the barycenter of the face, orthogonal to the face of the simplex, controlling the distance by curvature estimates. However, it introduces the possibility that the complex closes on itself and a new simplex overlaps an old.



[1]Eugene L. Allgower, Stefan Gnutzmann, An Algorithm for Piecewise Linear Approximation of Implicitly Defined Two-Dimensional Surfaces. SIAM Journal on Numerical Analysis, Volume 24, Number 2, 452--469, 1987.

 [2] E. L. Allgower, K. Georg, Simplicial and Continuation Methods for Approximations, Fixed Points and Solutions to Systems of Equations SIAM Review, Volume 22, 28--85, 1980.

 [3] Eugene L. Allgower, Phillip H. Schmidt, An Algorithm for Piecewise-Linear Approximation of an Implicitly Defined Manifold SIAM Journal on Numerical Analysis, Volume 22, Number 2, 322--346, April 1985.

[4] David P. Dobkin, Silvio V. F. Levy, William P. Thurston and Allan R. Wilks, Contour Tracing by Piecewise Linear Approximations. ACM Transactions on Graphics, 9(4) 389-423, 1990.

Predictor-Corrector Continuation Methods

A second class of continuation methods, predictor-corrector continuation, uses a mapping from R^k onto M at a point u_m(s), u_m(0)=u_m, and  A_m(u_m) is the image of a neighborhood of the origin under this map. The Implicit Function Theorem says that a unique mapping exists provided that u_m is a point at which the Jacobian F_u is full rank (in this case n-k). The mapping is only unique in the sense that there is a unique surface through u_m, many different parameterizations are possible. One approach is to find an orthonormal basis for the null space of F_u, and an iterative scheme to project orthogonally to the tangent space onto M. This is called the pseudo-arclength parameterization, or the exponential map.

For k=1, predictor-corrector continuation is well defined. It has not been easy to extend it to higher dimensions, basically because there is noway to introduce a  reference decomposition which ensures compatibility between the simplices being merged.

In one dimension, M_m is represented as an open polygonal arc with vertices on M. That is, the image of a regular grid on R. The point u_m is one endpoint of the arc. The neighborhood A_m(u_m) is the short arc [u_m(-ds_m),u_m,u_m(ds_m)]. What makes the continuation work when k=1 is that one interval of this arc ([u_m(-ds_m),u_m] if the tangent is always outward pointing), is always entirely interior to M_m. The merge then simply drops that interval, and unless the neighborhood contains the starting point u_0, the second interval is appended to the polygonal arc. If the starting point does lie in the interval, M has the same topology as the circle, and the second interval [u_m,u_m(ds_m)] is shortened to [u_m,u_0], and the continuation terminates.

One dimensional continuation as adding an interval (line segment)

This is approximating A_m(u_m)=u_m([-ds_m,ds_m]), and M_m(s)=u_i(s-s_i), where s is in the interval [s_i,s_{i+1}], s_i=s_{i-1}+ds_i, s_0=0.


[5] Willy J. F. Govaerts, Numerical Methods for Bifurcations of Dynamical Equilibria, SIAM. 2000.

[6] H. B. Keller, Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems, in Applications of Bifurcation Theory, P. Rabinowitz ed., Academic Press, 1977.

[7] W.C. Rheinboldt and J.V. Burkardt, A Locally Parameterized Continuation Process, ACM Transactions on Mathematical Software, Volume 9, 236--246, 1983.

[8] Morgan, Alexander, Solving polynomial systems using continuation for engineering and scientific problems, Prentice-Hall, Englewood Cliffs, N.J. 1987

[9] C. B. Garcia and W. I. Zangwill, Pathways to Solutions, Fixed Points and Equilibria. Prentice-Hall, 1981.

[10] E. Doedel, Nonlinear Numerics. International Journal of Bifurcation and Chaos, 7(9):2127-2143, 1997.

[11] R. Seydel, Nonlinear Computation. International Journal of Bifurcation and Chaos, 7(9):2105-2126, 1997.

Predictor-Corrector Algorithms in Higher Dimensions

Rheinbolt's First Algorithm

The first attempt to produce a multidimensional predictor-corrector continuation method, due to Rheinbolt, mimics the case of k=1, using a reference decomposition of R^k, and computing points on M corresponding to each vertex. The point u_m is a vertex on the boundary of M_m, and the neighborhood A_m(u_m) is the image of the cells in the reference decomposition which contain the boundary vertex. Some of these will be part of M_m, and some will need to be assigned to point on M. For k=1 the reference decomposition of the real line is a uniform mesh, the neighborhood is the pair of intervals [s_{m-1},s_m] and [s_m,s_{m+1}], on either side of the endpoint, and the point s_{m+1} is assigned to u_{m+1}.

The trick to getting this to work for k>1, is to ensure that u_m(s_i-s_m)=u_i for the points which have been given values. The unassigned vertices are found by evaluating the mapping. When k=1 this is that u_m(-ds_m)=u_{m-1}. Obviously, there is always a way to choose the parameterization of u_m(s) for which this is true. For higher dimensions it is sometimes possible to choose such a parameterization, and that is what Rheinbolt's first algorithm does.

Rheinbolt's algorithm as adding a teighborhood of triangles to a regular grid

The algorithm essentially finds a global Euclidean parameterization for M, which of course is not possible unless M is isomorphic to Euclidean k-space. Also, some choices for the local parameterizations may lead to a situation where a parameterization cannot be found, even though a global parameterization does exist. Finally, nothing has been said about the analogous case to u_0 lying in the neighborhood A_m(u_m). I call this global self-intersection, and it occurs when some point on the boundary of M_m lies in A_m(u_m). It indicates that the topology of the boundary of M_m is different from the boundary of M_{m+1}, and if not detected will cause the continuation to repeatedly recover M instead of terminating.


[14] W.C. Rheinbolt, On a Moving Frame Algorithm and the Triangulation of Equilibirum Manifolds, In T. Kuper, R. Seydel, and H. Troger eds. "ISNM79: Bifurcation: Analysis, Algorithms, Applications", pages 256-267. Birkhauser, 1987.

[15] W.C. Rheinbolt, On the Computation of Multi-Dimensional Solution Manifolds of Parameterized Eqautions. Numerishe Mathematik, 53, 1988, pages 165-181.

Rheinbolt and Brodzik's Algorithm

The problem with choosing the local parameterizations in Rheinbolt's first algorithm can be avoided if the reference decomposition is adapted as the continuation progresses, instead of being fixed. This is what is done in the two dimensional version of Rheinbolt and Brodzik's algorithm, and Brodzik's extension to higher dimensions. The neighborhood A_m(u_m) is
replaced with the image of a standard decomposition of the interior of a sphere in k dimensions. The manifold M_m is required to be local Delaunay, that is, if the simplices near u_i (measured on the complex) are projected into the tangent space to M at u_i, a circumsphere around any simplex contains no vertices of M_m other than the vertices of the simplex.

Brodzik's algorithm as adding a triangulation of a sphere

The simplices in A_m(u_m) (from sphere) and those in M_m won't usually be compatible. Instead of choosing the parameterization so that they are, the part of M_m near u_m is projected onto the tangent space, giving a set of simplices, which must be merged with those from the interior of the sphere. Only the vertices from the sphere are kept, and those which lie inside the circumsphere of any of the simplices from M_m are also discarded. This ensures that a Delaunay triangulation of the vertices from M_m and those remaining from  the sphere is compatible with M_m, and that when merged, M_{m+1} is locally Delaunay.

Brodzik, part 2. Merging the sphere with the existing triangulation

Finally, the remaining vertices of the sphere are projected using u_m(s), and included in M_{m+1}.

This removes the problem of choosing the parameterization, but as stated does not handle global self-intersection. If all of the simplices of M_m which lie near u_m in the embedding space R^n are projected onto the tangent space at u_m, this could be taken care of as well.


[16] M. L. Brodzik and W.C. Rheinbolt, On the Simplicial Approximation of Implicitly Defined Two-Dimensional Manifolds. Computers and Mathematics with Applications, 28(9): 9-21, 1994.

[17] M. L. Brodzik, The Computation of Simplicial Approximations of Implicitly Defined p-Manifolds. Computers and Mathematics with Applications, 36(6):93-113, 1998.

Melville and Mackey

Melville and Mackey proposed an algorithm for k=2 which avoids the problem of the merge by using a boundary representation for M_m and A_m(u_m). They are represented as the interior of a polygonal curve with vertices on M, instead of as a triangulation of the interior. The merge then becomes a question of finding the polygonal curve which is the boundary of the interior of two curves, which can be done easily by walking around the curves checking for crossings, and dividing intervals which cross into two non-crossing intervals, then discarding whole intervals. This essentially reproduces the advantage gained from the representation as polygonal arc from the k=1 case to k=2.

Melville amd <ackey's algorithm as merging boundary paths.


[18] R. Melville and D. S. Mackey, New Algorithm for Two-Dimensional Numerical Continuation. Computers and Mathematics with Applications, 30(1):31-46, 1995.


Here at Watson we are developing predictor-corrector methods for computing higher dimensional implicitly defined manifolds, using a different generalization of the k=1 algorithm. The pair of intervals [-ds_m,0],[0,ds_m] might be viewed as a "triangulation" of part of the real line, which leads to the algorithms described above, or as a "ball", |s|<ds_m. This second generalization turns out to lead to a much simpler algorithm, reducing the algorithm to the operation of updating the vertices of a polyhedron in k space as half spaces are removed from it.

 An older version uses an ellipse in the tangent plane to compute surfaces. There is a research report available which describes it. The new version is available as multifario, and  we have some examples.

Basically, if the neighborhoods are A_m(u_m)=u_m(B_R_i(0)), (which is u_m(s), |s|<R_i) and M_m is the union of the neighborhoods, the merge operation is trivial, and the difficulty is to find a point on the boundary of M_m. If M were flat this means finding the boundary of a set of spherical balls of various radii. The boundary is made up of pieces of the spheres, precisely those that lie in none of the others, that is

d M_m = U (dB_R_i(s_i) - U dB_R_i(s_i)^B_R_j(s_j))

Now, the repeated subtractions can be written as the intersection of pairwise subtractions:

dB_R_i(s_i) - U dB_R_i(s_i)^B_R_j(s_j) = ^ (dB_R_i(s_i) -dB_R_i(s_i)^B_R_j(s_j))

and for spheres this is particularly simple. The part of a sphere that lies outside of another sphere is the part in a halfspace perpendicular to the lie connecting the centers of the spheres. This is true in any dimension (symmetry!) and means that the expression above is equivalent to the intersection of the surface of a sphere and a convex (possibly infinite) polyhedron. If that polyhedron is clipped to a cube containing the sphere the expression is unchanged (the cube contains the sphere). And that means that the boundary of the union of spherical balls is the union of the intersection of the boundary of each ball with a finite, convex polyhedron which is the intersection of a cube with a set of half spaces defined by the overapping balls. (That's a compcated statement, but if you think on it, it's obvious.)


  [19] Henderson, Michael E., Multiple Parameter Continuation: Computing Implicitly Defined k-manifolds, IJBC v12(3), pages 451-76.


The examples above were nice in that the solution curve cut each triangle at exactly two places, and there was a unique tangent at each point on the curve. All sorts of interesting things happen when these don't hold. One curve will typically split into several curves at such points, and so they are called bifurcation points.

 The basic idea of bifurcation theory is that locally the solution curve can be made to coincide with solutions of an algebraic system by non-singular transformations. These algebraic models then become a classification for the possible behavior of a system near a singular point.

 The classification begins with the Jacobian of F and the left and right null spaces of the Jacobian. If the Jacobian at a point u is full rank (the right null space is dimension k and the left null space is dimension 0), then the Implicit Function Theorem states that the solution is the image of a linear k-space (there is a u(s) such that F(u(s))=0).

 If the Jacobian is not full rank, then a decomposition is used to project the equations onto the left null space and the unknowns onto the right null space. This is called the Lyapunov-Schmidt decomposition. If Phi={Phi_1,....,Phi_m} is an orthonormal basis for the right null space, and Psi={Psi_1,...,Psi_n} is an orthonormal basis for the left null space, then the decomposition is

A Lyapunov-Schmidt decomposition using the left and right null spaces.

Here Phi projects into the right null space of F_u and Psi projects into the left null space of F_u. Therefore, with

Decomposition of independant variabls u=(xi,a)

F(u)=0 is equivalent to the pair of equations

Decomposition of F=(G(xi,a),H(xi,a))

The first of these has Jacobian wrt to a (a little abuse of the notation here)

G_a=(I-Psi Psi^T)F_u (I-Phi Phi^T)

which is nonsinglar (the inverse is the Pseudoinverse of F_u). The Implicit Function Theorem can be used to solve this for a as a function of xi, leaving the bifuration equations

B(xi)=H(xi,a(xi))=Phi^T F(u_0+Phi xu +a(xi))=0

The bifurcation equations are still nonlinear, but a much smaller system than the original system. In some sense the decomposition is concentrating the nonlinearity into the bifurcation equations and removing the linearity. A Taylor series for the bifurcation equation starts:

B(xi) sim 1/2 |si^T F_uu(u_0) Phi Phi xi xi + ...

Where the first few derivatives of a are:

a(0)=0, a_xi(0)=0, ...

A classification of solutions of the bifurcation equations can be performed, by asking how many of the higher order terms can be removed without changing the solution. With

A^i_{jk}=\Psi_i^T F_{uu} phi_j phi_k, B^i_{jkl}=...

This results in the following first few cases:

Name n m Model Requirements Diagram
Quadratic bifurcation 1 2 xi_0^2-xi_1^2=0 A^0_{01}A^0_{01}-A^0_{00}A^2_{11}>0 Crossing transcritical solution curves
Isola center 1 2 xi_0^2-xi_2 A^0_{01}A^0_{01}-A^0_{00}A^0_{11}<0 The single point solution Isola center manifold

The classification procedes first in terms of the dimensions of the left and right null spaces of the Jacobian, then in terms of the solutions of the algebraic bifurcation equations. For low orders and a single algebraic equation this is feasible, classifying the solutions of higher order polynomials in terms of their coefficients is more problematic. Classifying the solutions of an algebraic system is done only in special cases. One question of computational bifurcation theory is how to write systems of equations for singular points that are themselves non-singular. This has been done for many types of singular point, and is used to study systems which involve more than a single parameter.

Perturbed Bifurcation Theory

Some of these are stable to general perturbations, others are not. The analysis basically looks for a higher dimensional system which is stable to perturbation, and which has the unstable configuration as a section. This is called unfolding a singularity. For example, the Cubic bifurcation is unstable to general perturbation (though stable to ``symmetric'' perturbations), and is a section of the Cusp catastrophe. This leads to another classification, called Catastrophe Theory.

Name Model
Fold catastrophe x^2 + s = 0
Cusp catastrophe x^3 + s x + t = 0
Swallowtail catastrophe x^4 + s x^2 + t x + u = 0
Butterfly catastrophe x^5 + s x^3 + t x^2 + u x + v = 0

Symmetry and Bifurcation Theory

The perturbations described above are general. One school of thought is that everything in nature is perturbed, and therefore we shouldn't observe phenomena which are not stable. This may be true, but some perturbations in nature preserve some symmetry of the problem, and therefore we also need to ask what happens in the presence of symmetries.

Detecting Bifurcation and Branch Switching

Bifurcation of Dynamical Systems

For parameterized dynamical systems, it is solutions of equations like -

u_t=F(u,lambda,t) or u_{i+1) = F(u_i,lambda)

that is of interest. Steady states of these are the case described above, but other interesting things can happen in the larger context.

 Steady states can bifurcate to periodic motions (Hopf bifurcation), periodic motions may period double (or other multiples), infinite period motions (homoclinic and heteroclinic trajectories) can bifurcate from steady state motions (Takens-Bogdanov) and whole families of periodic and homoclinic and heteroclinic motions can bifurcate at a Shilnikov bifurcation.

 When computing steady states and bifurcations from steady states to other types of solutions, the analysis is as before, but other tests (beyond checking for a real eigenvalue crossing zero) must be used to detect bifurcations. For Hopf bifurcation this is a complex pair of eigenvalues crossing the imaginary axis. For period multiplying in discrete maps it is eigenvalues crossing the unit circle.

 To compute solutions that are not steady states a system is formulated for which they are fixed points. For discrete maps this is folding the function, i.e. a doubled point satisfies --

u_1=F(u_0,lambda), u_0=F(u_1,lambda)

So we look for steady states of


For periodic motions (Poincar\'e Continuation) in the continuous systems the curve of u(t) is discretized on [0,1] and a phase constraint is introduced --

u_t = T F(u(t),lambda), u(0)=u(1), phase constraint

where u_0(t) is a reference periodic motion, usually the solution at a nearby parameter value.

 Other structures have been formulated for computation, including homoclinic and heteroclinic connections, and invariant tori.

Computing Manifolds of Singular Points

[Research home page]

 [ IBM home page | Order | Privacy | ContactIBM | Legal ]