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:

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  of a point  on the manifold is merged with the previously computed part of the manifold , which produces . A new point  is then selected from , and the process is repeated until the part of M inside some finite region is covered. The continuation starts with a single point , provided by the user, and as long as has some part which is exterior to the process will eventually terminate. The natural choice is therefore that  be on or near the boundary of .

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 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 . This is what simplicial continuation does. It represents  as the intersection of with an n-k dimensional simplex . The neighborhood  is  then the intersection of  and a set of n dimensional simplices with  in its interior. Of course, the intersection need not be done explicitly, it is enough just to have  be the n dimensional simplices. If   is represented in the same way, then if the simplices in  and  are compatible, (the intersection of two simplcies is either empty or another simplex), then the merge operation is just adding the simplices from  which are not already present. To ensure this is the case,  and  are subsets of a reference simplicial decomposition of .

The next ''point'', , is any n-k dimensional face on the boundary of  which crosses .
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  for crossing M may require quite a bit of work if n is large (only those exterior to 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 assigns a value in to each vertex of a (n-k)-dimensional simplex in . So for example, a scalar to each endpoint of an edge in the plane-

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 , 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 --

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,

,

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(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

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  to  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  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.

Triangulations

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  (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.

References:

[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  at a point u_m(s), u_m(0)=u_m, and   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 , 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 . 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,  is represented as an open polygonal arc with vertices on . That is, the image of a regular grid on R. The point is one endpoint of the arc. The neighborhood  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 . The merge then simply drops that interval, and unless the neighborhood contains the starting point , the second interval is appended to the polygonal arc. If the starting point does lie in the interval,  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.

This is approximating =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.

References:

[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 , and the neighborhood  is the image of the cells in the reference decomposition which contain the boundary vertex. Some of these will be part of , 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 .

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.

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 . I call this global self-intersection, and it occurs when some point on the boundary of lies in . It indicates that the topology of the boundary of  is different from the boundary of M_{m+1}, and if not detected will cause the continuation to repeatedly recover M instead of terminating.
 

References

[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 is
replaced with the image of a standard decomposition of the interior of a sphere in k dimensions. The manifold 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  other than the vertices of the simplex.

The simplices in  (from sphere) and those in  won't usually be compatible. Instead of choosing the parameterization so that they are, the part of  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  and those remaining from  the sphere is compatible with , and that when merged, M_{m+1} is locally Delaunay.

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  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.

References

[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  and . 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.

References

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

Henderson

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  is the union of the neighborhoods, the merge operation is trivial, and the difficulty is to find a point on the boundary of . 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

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

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.
 

Bifurcation

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  is an orthonormal basis for the right null space, and  is an orthonormal basis for the left null space, then the decomposition is

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

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

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

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

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:

Where the first few derivatives of a are:

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

This results in the following first few cases:

Name	n	m	Model	Requirements	Diagram
Quadratic bifurcation	1	2
Isola center	1	2

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
Cusp catastrophe
Swallowtail catastrophe
Butterfly catastrophe

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 -

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 --

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 --

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 ]