# Numerical continuation

﻿
Numerical continuation

Numerical continuation is a method of computing approximate solutions of a system of parameterized nonlinear equations, $F(\mathbf u,\lambda) = 0$

The parameter λ is usually a real scalar, and the solution an n-vector. For a fixed parameter value λ,, $F(\ast,\lambda)$ maps Euclidean n-space into itself.

Often the original mapping F is from a Banach space into itself, and the Euclidean n-space is a finite dimensional approximation to the Banach space.

A steady state, or fixed point, of a parameterized family of flows or maps are of this form, and by discretizing trajectories of a flow or iterating a map, periodic orbits and heteroclinic orbits can also be posed as a solution of F = 0.

## Other forms

In some nonlinear systems, parameters are explicit. In others they are implicit, and the system of nonlinear equations is written $F(\mathbf u) = 0$

where $\mathbf u$ is an n-vector, and its image $F(\mathbf u)$ is an n-1 vector.

This formulation, without an explicit parameter space is not usually suitable for the formulations in the following sections, because they refer to parameterized autonomous nonlinear dynamical systems of the form: $\mathbf u' = F(\mathbf u,\lambda)$

However, in an algebraic system there is no distinction between unknowns $\mathbf u$ and the parameters.

## Periodic motions

A periodic motion is a closed curve in phase space. That is, for some period T. $\mathbf u' = F(\mathbf u,\lambda),\, \mathbf u(0) = \mathbf u(T)$

The textbook example of a periodic motion is the undamped pendulum.

If the phase space is periodic in one or more coordinates, say $\mathbf u(t) = \mathbf u(t) + \Omega$, with Ω a vector, then there is a second kind of periodic motions defined by $\mathbf u' = \mathbf F(\mathbf u,\lambda),\, \mathbf u(0) = \mathbf u(T) + N.\Omega$

Here N is a diagonal matrix of integers that serves as an index of these periodic motions of the second kind.  The first step in writing an implicit system for a periodic motion is to move the period T from the boundary conditions to the ODE: $\mathbf u' = T\mathbf F(\mathbf u,\lambda),\, \mathbf u(0)=\mathbf u(1) + N.\Omega$

The second step is to add an additional equation, a phase constraint, that can be thought of as determining the period. This is necessary because any solution of the above boundary value problem can be shifted in time by an arbitrary amount (time does not appear in the defining equations—the dynamical system is called autonomous).

There are several choices for the phase constraint. If $\mathbf u_0(t)$ is a known periodic orbit at a parameter value λ0 near λ, then, Poincaré used $<\mathbf u(0) - \mathbf u_0(0),\mathbf F(\mathbf u_0(0),\lambda_0)> = 0$

which states that $\mathbf u$ lies in a plane which is orthogonal to the tangent vector of the closed curve. This plane is called a Poincaré section. For a general problem a better phase constraint is an integral constraint introduced by Eusebius Doedel, which chooses the phase so that the distance between the known and unknown orbits is minimized: $\int_0^1 <\mathbf u(t) - \mathbf u_0(t),\mathbf F(\mathbf u_0(t),\lambda_0)>\, dt = 0$

## Homoclinic and heteroclinic motions  ## Definitions

### Solution component

A solution component $\Gamma(\mathbf u_0,\lambda_0)$ of the nonlinear system F is a set of points $(\mathbf u,\lambda)$ which satisfy $F(\mathbf u,\lambda)=0$ and are connected to the initial solution $(\mathbf u_0,\lambda_0)$ by a path of solutions $(\mathbf u(s),\lambda(s))$ for which $(\mathbf u(0),\lambda(0))=(\mathbf u_0,\lambda_0),\, (\mathbf u(1),\lambda(1)) = (\mathbf u,\lambda)$

and $F(\mathbf u(s),\lambda(s))=0$. This figure shows two solution components, one red and the other blue. Note that these two components may be connected outside the region of interest.

### Numerical continuation

A numerical continuation is an algorithm which takes as input a system of parametrized nonlinear equations and an initial solution $(\mathbf u_0,\lambda_0)$, $F(\mathbf u_0,\lambda_0)=0$, and produces a set of points on the solution component $\Gamma(\mathbf u_0,\lambda_0)$.

### Regular point

A regular point of F is a point $(\mathbf u,\lambda)$ at which the Jacobian of F is full rank (n).

Near a regular point the solution component is an isolated curve passing through the regular point (the implicit function theorem). In the figure above the point $(\mathbf u_0,\lambda_0)$ is a regular point.

### Singular point

A singular point of F is a point $(\mathbf u,\lambda)$ at which the Jacobian of F is not full rank.

Near a singular point the solution component may not be an isolated curve passing through the regular point. The local structure is determined by higher derivatives of F. In the figure above the point where the two blue curves cross is a singular point.

In general solution components Γ are branched curves. The branch points are singular points. Finding the solution curves leaving a singular point is called branch switching, and uses techniques from bifurcation theory (singularity theory, catastrophe theory).

For finite dimensional systems (as defined above) the Lyapunov-Schmidt decomposition may be used to produce two systems to which the Implicit Function Theorem applies. The Lyapunov-Schmidt decomposition uses the restriction of the system to the complement of the null space of the Jacobian and the range of the Jacobian.

If the columns of the matrix Φ are an orthonormal basis for the null space of $J=\left[ \begin{array}{cc} F_x & F_{\lambda}\\ \end{array} \right]\,$ and the columns of the matrix Ψ are an orthonormal basis for the left null space of J, then the system F(x,λ) = 0 can be rewritten as $\left[ \begin{array}{l} (I-\Psi\Psi^T)F(x+\Phi\xi + \eta)\\ \Psi^T F(x+\Phi\xi + \eta)\\ \end{array} \right]=0,$ where η is in the complement of the null space of J $(\Phi^T\,\eta=0)$.

In the first equation, which is parametrized by the null space of the Jacobian (ξ), the Jacobian with respect to η is non-singular. So the implicit function theorem states that there is a mapping η(ξ) such that η(0) = 0 and (I − ΨΨT)F(x + Φξ + η(ξ)) = 0). The second equation (with η(ξ) substituted) is called the bifurcation equation (though it may be a system of equations).

The bifurcation equation has a Taylor expansion which lacks the constant and linear terms. By scaling the equations and the null space of the Jacobian of the original system a system can be found with non-singular Jacobian. The constant term in the Taylor series of the scaled bifurcation equation is called the algebraic bifurcation equation, and the implicit function theorem applied the bifurcation equations states that for each isolated solution of the algebraic bifurcation equation there is a branch of solutions of the original problem which passes through the singular point.

Another type of singular point is a turning point bifurcation, or saddle-node bifurcation, where the direction of the parameter λ reverses as the curve is followed. The red curve in the figure above illustrates a turning point.

## Particular algorithms

### Natural parameter continuation

Most methods of solution of nonlinear systems of equations are iterative methods. For a particular parameter value λ0 a mapping is repeatedly applied to an initial guess $\mathbf u_0$. If the method converges, and is consistent, then in the limit the iteration approaches a solution of $F(\mathbf u,\lambda_0)=0$.

Natural parameter continuation is a very simple adaptation of the iterative solver to a parametrized problem. The solution at one value of λ is used as the initial guess for the solution at λ + Δλ. With Δλ sufficiently small the iteration applied to the initial guess should converge. One advantage of natural parameter continuation is that is uses the solution method for the problem as a black box. All that is required is that an initial solution can be given (some solvers used to always start at a fixed initial guess). There has been a lot of work in the area of large scale continuation on applying more sophisticated algorithms to black box solvers. (see e.g. LOCA).

However, natural parameter continuation fails at turning points, where the branch of solutions turns round. So for problems with turning points, a more sophisticated method such as pseudo-arclength continuation must be used (see below).

### Simplicial or piecewise linear continuation

Simplicial Continuation, or Piecewise Linear Continuation (Allgower and Georg) is based on three basic results.

The first is

 If F(x) maps IR^n into IR^(n-1), there is a unique linear interpolant on an (n-1)-dimensional simplex which agrees with the function values at the vertices of the simplex.

The second result is:

 An (n-1)-dimensional simplex can be tested to determine if the unique linear interpolant takes on the value 0 inside the simplex.

Please see the article on piecewise linear continuation for details.

With these two operations this continuation algorithm is easy to state (although of course an efficient implementation requires a more sophisticated approach. See [B1]). An initial simplex is assumed to be given, from a reference simplicial decomposition of IR^n. The initial simplex must have at least one face which contains a zero of the unique linear interpolant on that face. The other faces of the simplex are then tested, and typically there will be one additional face with an interior zero. The initial simplex is then replaced by the simplex which lies across either face containing zero, and the process is repeated. References: Allgower and Georg [B1] provides a crisp, clear description of the algotihm.

### Pseudo-arclength continuation

This method, which was proposed by H.B. Keller [A5] in the late 1970s, is based on the observation that the "ideal" parameterization of a curve is arclength. Pseudo-arclength is an approximation of the arclength in the tangent space of the curve. The resulting modified natural continuation method makes a step in pseudo-arclength (rather than λ). The iterative solver is required to find a point at the given pseudo-arclength, which requires appending an additional constraint (the pseudo-arclength constraint) to the n by n+1 Jacobian. It produces a square Jacobian, and if the stepsize is sufficiently small the modified Jacobian is full rank. The algorithm is a predictor-corrector method. The prediction step finds the point (in IR^(n+1) ) which is a step Δs along the tangent vector at the current pointer. The corrector is usually Newton's method, or some variant, to solve the nonlinear system $\begin{array}{l} F(u,\lambda)=0\\ \dot u^*_0(u-u_0)+\dot \lambda_0 (\lambda-\lambda_0) = \Delta s\\ \end{array}\,$

where $(\dot u_0,\dot\lambda_0)\,$ is the tangent vector at $(u_0,\lambda_0)\,$. The Jacobian of this system is the bordered matrix $\left[ \begin{array}{cc} F_u & F_{\lambda}\\ \dot u^* & \dot \lambda\\ \end{array} \right]\,$

At regular points, where the unmodified Jacobian is full rank, the tangent vector spans the null space of the top row of this new Jacobian. Appending the tangent vector as the last row can be seen as determining the coefficient of the null vector in the general solution of the Newton system (particular solution plus an arbitrary multiple of the null vector).

### Gauss–Newton continuation

This method is a variant of pseudo-arclength continuation. Instead of using the tangent at the initial point in the arclength constraint, the tangent at the current solution is used. This is equivalent to using the pseudo-inverse of the Jacobian in Newton's method, and allows longer steps to be made. [B17]

## Continuation in more than one parameter

The parameter λ in the algorithms described above is a real scalar. Most physical and design problems generally have many more than one parameter. Higher dimensional continuation refers to the case when λ is a k-vector.

The same terminology applies. A regular solution is a solution at which the Jacobian is full rank (n). A singular solution is a solution at which the Jacobian is less than full rank.

A regular solution lies on a k-dimensional surface, which can be parameterized by a point in the tangent space (the null space of the Jacobian). This is again a straightforward application of the Implicit Function Theorem.

## Applications of numerical continuation techniques

Numerical continuation techniques have found a great degree of acceptance in the study of chaotic dynamical systems and various other systems which belong to the realm of catastrophe theory. The reason for such usage stems from the fact that various non-linear dynamical systems behave in a deterministic and predictable manner within a range of parameters which are included in the equations of the system. However, for a certain parameter value the system starts behaving chaotically and hence it become necessary to follow the parameter in order to be able to decipher the occurrences of when the system starts being non-predictable, and what exactly (theoretically) makes the system become unstable.

Analysis of parameter continuation can lead to more insights about stable/critical point bifurcations. Study of saddle-node, transcritical, pitch-fork, period doubling, Hopf, secondary Hopf (Neimark) bifurcations of stable solutions allows for a theoretical discussion of the circumstances and occurrences which arise at the critical points. Parameter continuation also gives a more dependable system to analyze a dynamical system as it is more stable than more interactive, time-stepped numerical solutions. Especially in cases where the dynamical system is prone to blow-up at certain parameter values (or combination of values for multiple parameters).

It is extremely insightful as to the presence of stable solutions (attracting or repelling) in the study of Nonlinear Partial Differential Equations where times stepping in the form of the Crank Nicolson algorithm is extremely time consuming as well as unstable in cases of nonlinear growth of the dependent variables in the system. The study of turbulence is another field where the Numerical Continuation techniques have been used to study the advent of turbulence in a system starting at low Reynold's numbers. Also, a research using these techniques has provided the possibility of finding stable manifolds and bifurcations to invariant-tori in the case of the restricted three body problem in Newtonian gravity and have also given interesting and deep insights into the behaviour of systems such as the Lorenz equations.

## Software

(Under Construction) See also The SIAM Activity Group on Dynamical Systems' list http://www.dynamicalsystems.org/sw/sw/

## Examples

This problem, of finding the points which F maps into the origin appears in computer graphics as the problems of drawing contour maps (n=2), or isosurface(n=3). The contour with value h is the set of all solution components of F-h=0

Wikimedia Foundation. 2010.

### Look at other dictionaries:

• Continuation War — Part of Eastern Front of World War II Finnish StuG III Ausf. G assault guns on parade …   Wikipedia

• Numerical sign problem — The numerical sign problem refers to the difficulty of numerically evaluating the integral of a highly oscillatory function of a large number of variables. Numerical methods fail because of the near cancellation of the positive and negative… …   Wikipedia

• Piecewise linear continuation — implicial ContinuationSimplicial Continuation, or Piecewise Linear Continuation (Allgower and Georg  ,  ) is a one parameter continuation method which is well suited to small to medium embedding spaces. The algorithm has been generalized to …   Wikipedia

• List of numerical analysis topics — This is a list of numerical analysis topics, by Wikipedia page. Contents 1 General 2 Error 3 Elementary and special functions 4 Numerical linear algebra …   Wikipedia

• List of mathematics articles (N) — NOTOC N N body problem N category N category number N connected space N dimensional sequential move puzzles N dimensional space N huge cardinal N jet N Mahlo cardinal N monoid N player game N set N skeleton N sphere N! conjecture Nabla symbol… …   Wikipedia

• Predictor-corrector method — In mathematics, particularly numerical analysis, a predictor corrector method is an algorithm that proceeds in two steps. First, the prediction step calculates a rough approximation of the desired quantity. Second, the corrector step refines the… …   Wikipedia

• Fortsetzungsverfahren — Homotopie Verfahren (auch als Homotopiemethode, Fortsetzungs oder Einbettungsverfahren bezeichnet) sind Berechnungsmethoden in der numerischen Mathematik zur Bestimmung von Lösungen nichtlinearer Gleichungssysteme. Ziel ist es dabei den… …   Deutsch Wikipedia

• Homotopieverfahren — Homotopie Verfahren (auch als Homotopiemethode, Fortsetzungs oder Einbettungsverfahren bezeichnet) sind Berechnungsmethoden in der numerischen Mathematik zur Bestimmung von Lösungen nichtlinearer Gleichungssysteme. Ziel ist es dabei den… …   Deutsch Wikipedia

• Speed Demon (comics) — Superherobox| caption=Speed Demon battles Spider Man on the cover of Amazing Spider Man #222 (Nov. 1981). Art by Bob Hall. character name=Speed Demon real name=James Sanders publisher=Marvel Comics debut= (as the Whizzer) Avengers #70 (Nov. 1969) …   Wikipedia

• KABBALAH — This entry is arranged according to the following outline: introduction general notes terms used for kabbalah the historical development of the kabbalah the early beginnings of mysticism and esotericism apocalyptic esotericism and merkabah… …   Encyclopedia of Judaism