New: My continuation code "Multifario" is now available under an Open Source license from COIN-OR.
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 consist of ``regular'' pieces, which are joined at ``singular'' solutions (which are also solution manifolds, but of a system with different "n" and "k-1"). The regular pieces are curves when k=1, surfaces when k=2 and k-manifolds in general. These systems 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.
Hinke Osinga has a page of links to software
for the study of dynamical systems that includes continuation
software. DSWEB
has lots of good information including a list of tutorials on dynamical
systems. R. Seydel has a set of tutorials and lots of examples on Bifurcation and Continuation Methods.
The basic idea of a continuation method is very simple: compute a
piece
of the solution manifold near one solution, then select another
solution
from this set and repeat the process. As long as the new piece covers
some
new part of the solution manifold the computation progresses.So the
basic
issues are
There are two ways to do the first task, and these divide continuation methods into two types: 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.
The different methods make different choices for how the neighborhoods are computed, and how the point and manifold are represented.
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-
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 --
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,
,
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.
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.
[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.
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.
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.
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 Rheinboldt'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 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.
[15] W.C. Rheinboldt, On the Computation of Multi-Dimensional Solution Manifolds of Parameterized Eqautions. Numerishe Mathematik, 53, 1988, pages 165-181.
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.
[17] M. L. Brodzik, The Computation of Simplicial Approximations of Implicitly Defined p-Manifolds. Computers and Mathematics with Applications, 36(6):93-113, 1998.
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.
Basically, if the neighboroods 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
Now, the repeated subtractions can be written as the intersection of pairwise subtractions:
and for spheres this is particularly simple.
This result says that to find the part of the sphere that lies
on the boundary we must take the intersection of the intersection of
the
sphere with these half-spaces. In fact we can first intersect the
sphere
with a cube containing it, and so the part of the sphere that lies on
the
boundary is the part in the convex polyhedron formed by removing
halfspaces
from a cube containing the sphere, one halfspace for each overlapping
sphere.
 |
 |
The algorithm is this:
Now, this is really an algorithm for computing a "flat" space. For
curved
manifolds the balls are in the tangent spaces at the centers. To do the
subtraction, we project the neighboring balls into the same tangent
space,
and approximate them by balls -- in the paper it is shown that this is
valid as long as the distance between the tangent plane and the
manifold
is small inside the ball.
 |
 |
 |
That's it. We have reduced the problem of finding of a boundary point to updating the vertices of a set of polyhedra as half spaces are removed from a cube. This is a straight forward operation. If a list of the polyhedra which might have an exterior vertex is kept, and a binary tree of hierarchical bounding boxes for the balls (to make it easier to find overlapping balls), we get an N ln(N) continuation algorithm.
I want to point out that there is a direct analogy to what might be called "the advancing front approach" to this problem. If you think of Huygen's pronciple from optics and how it advances a wavefront:
This approach can also be used to construct an adaptive simplicial continuation method. If we choose the put the balls in n-space, instead of k-space, the representation of the manifold is OK (i.e. we can have represent the manifold as the union of the intersection of each polyhedron with the manifold), but the selection of the "next point" has to be modified. The n-k dimensional faces of the polyhedra will intersect the manifold at a single point, and the techniques desceribed above under "simplicial continuation" can be used to determine if such a face does contain an intersection point. So instead of projecting an exterior vertex we test each n-k dimensional face of the original polyhedron for transversality, and create a new ball at the intersection of a transverse face.
Here's a sketch of how I look at the simplicial version. The blue curve is the solution manifold, and the red the piecewise linear approximation. The black polyhedra are the "simplicies" (yeah, yeah, I know they're not simplicies). A black square marks vertices with F>0, and black circles F<0.
[1] Henderson, M.E,."Multiple Parameter Continuation: Computing Implicitly Defined k-manifolds", IJBC v12(3), pages 451-76
[2] Michael E. Henderson & Sébastien Neukirch, "Classification of the spatial equilibria of the clamped elastica: numerical continuation of the solution set", submitted to the IJBC. preprint
[3] Henderson, M.E., "Computing Implicitly Defined Surfaces: Two Parameter Continuation", IBM Research Report RC18777. preprint
and at which the Jacobian of F at u is full rank (n-k) are called regular points. In the introduction I asserted the solution set of the system is a k dimensional manifolds if all solutions are regular points. This is just a way of stating the Implicit Function Theorem (IFT). Another consequence of the IFT is that regular points cannot be on the boundary of the solution set (a full neighborhood of a regular solution exists). So regular ``branches'' cannot simply end, they must approach infinity, close upon themselves, or be bounded by a set of singular solutions (or a combination of these). This leads to the picture of the solution set as a number of pieces of k-dimensional manifold which are ``glued'' together at lower dimensional singular sets.
These bifurcations lead to solutions of the same system. Bifurcation theory also deals with bifurcation to solutions of equations for which F(u)=0 is a special case. For example, F(u)=0are steady states of the dynamical system u'=F(u). At a Hopf bifurcation a branch of periodic motions bifurcates from a branch of steady states. This type of bifurcation is much more interesting than the ones described below, are several good books are available on various aspects of it, including
Golubitsky and Schaeffer, Singularities and Groups in Bifurcation Theory, volume 1. Springer-Verlag, 1985.
Golubitsky, Stewart and Schaeffer, Singularities and Groups in Bifurcation Theory, volume2. Springer-Verlag, 1988.
Guckenheimer and Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer-Verlag, 1983.
Iooss and Joseph, Elementary Stability and Bifurcation Theory. Springer-Verlag 1980.
Kuznetsov, Elements of Applied Bifurcation Theory. Springer-Verlag, 1995.
Wiggins, Global Bifurcation and Chaos, Springer-Verlag, 1988.
The null space and range space of the Jacobian F_u define subspaces of n-space and (n-k)-space (the domain and range of F). One of the basic ideas of linear algebra is that there is an invertible mapping from the space orthogonal to the null space onto the range space. (The inverse of the mapping is the Pseudo-inverse of the Jacobian). This gives a natural splitting of both the domain and range of F, called the Lyapunov-Schmidt decomposition. If u* is a singular point, and F is a basis for the null space of the Jacobian at u*:
,
and 
,
,
and 
with 
.
.
.All this means that :
Now the inverse of the Jacobian of the second with respect to a is the Pseudoinverse of the Jacobian of F, so the Implicit Function Theorem says that there is a unique function
The coefficients in a Taylor series expansion for a can be found by repeated differentiation:
Notice that each derivative of a satisfies the same system with different right hand sides. In principle we can continue this as long as we have time and energy, but the expressions for the right hand side gets longer fast.
Now we take the bifurcation equations, and start differentiating:
Here F is a vector, F_u is a matrix, and F_uu and so on are tensors. When F is a scalar the first term in the bifurcation equation is usually written as a symmetric matrix with the same dimension as the null space of the Jacobian.
OK, so what happens at a singular point? Well, we need to solve the bifurcation equations. In general that's very difficult, but several cases have been dealt with explicitly, especially when the left null space is dimension one, and if symmetries exists the equations can be simplified.
Simple Bifurcation
Sometimes singular points with a single left null vector are called simple singular pointsand the analysis is somewhat simpler than the higher dimensional case. The second derivative can be written as a symmetric matrix:
=QvThen
Quadratic
Simple Bifurcation
A modified version of the Implicit Function can be used to find solutions near the origin. Basically we introduce a new variable, e, replace h by es, scale the bifurcation equation by 1/e^2, and introduce another equation |s|=1. (We've introduced one unknown and one equation.)
This
brings the
first term back down into the constant position, and if we look for a
solution
at e=0,
we find that s(0) must satisfy the Algebraic Bifurcation
Equations:
Let's look at a solution s of this and
consider the function B(e,s)=0 as a
function of e.
The IFT says then that if s*Q*B_xxx(0)QQss!=0, there is
a curve s(e) so that B(e,s(e))=0 in a
neigborhood of e=0,
and s(0)=s.
These are curves
on the solution manifold that pass through the singular point. If
the right null space is two
or three dimensional the cases are
 ,
Model |
 |
|
2,  |
h(e)=(0,0)
Isolated point |
 |
2,  |
h(e)=
 Crossing Curves |
 |
2,  |
h(e)=
 Isolated Curve |
 |
3,  |
h(e)= (0,0)
Isolated Point |
 |
3,  |
h(e,q)=
 Cone |
 |
3,  |
h(e,q)=
 Crossing Surfaces |
 |
3,  |
h(e,q)=
 Isolated Line |
 |
and so on. Each solution then is u=u*+FQh+a(Qh)
[ IBM home page | Order | Privacy | ContactIBM | Legal ]