Numerical continuation

Last updated

Numerical continuation is a method of computing approximate solutions of a system of parameterized nonlinear equations,

Contents

[1]

The parameter is usually a real scalar, and the solution an n-vector. For a fixed parameter value, maps Euclidean n-space into itself.

Often the original mapping is from a Banach space into itself, and the Euclidean n-space is a finite-dimensional 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 .

Other forms

In some nonlinear systems, parameters are explicit. In others they are implicit, and the system of nonlinear equations is written

where is an n-vector, and its image 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:

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

Periodic motions

A periodic motion is a closed curve in phase space. That is, for some period,

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

If the phase space is periodic in one or more coordinates, say , with a vector [ clarification needed ] , then there is a second kind of periodic motions defined by

for every integer .

PeriodicMotion.gif PeriodicOrbitTime.gif

The first step in writing an implicit system for a periodic motion is to move the period from the boundary conditions to the ODE:

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 is a known periodic orbit at a parameter value near , then, Poincaré used

which states that lies in a plane which is orthogonal to the tangent vector of the closed curve. This plane is called a Poincaré section .

PoincareSection.gif

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:

Homoclinic and heteroclinic motions

HomoclinicOrbit.gif HomoclinicOrbitTime.gif

Definitions

Solution component

Two solution components, one red and the other blue. Note that these two components may be connected outside the region of interest. Numerical Continuation components.png
Two solution components, one red and the other blue. Note that these two components may be connected outside the region of interest.

A solution component of the nonlinear system is a set of points which satisfy and are connected to the initial solution by a path of solutions for which and .

Numerical continuation

A numerical continuation is an algorithm which takes as input a system of parametrized nonlinear equations and an initial solution , , and produces a set of points on the solution component .

Regular point

A regular point of is a point at which the Jacobian of is full rank .

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 is a regular point.

Singular point

A singular point of is a point 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 . 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

and the columns of the matrix are an orthonormal basis for the left null space of , then the system can be rewritten as

where is in the complement of the null space of .

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 and . 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 a mapping is repeatedly applied to an initial guess . If the method converges, and is consistent, then in the limit the iteration approaches a solution of .

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.

NaturalParameter.gif

One advantage of natural parameter continuation is that it uses the solution method for the problem as a black box. All that is required is that an initial solution 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 . 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.

Simplicial.gif

References: Allgower and Georg [B1] provides a crisp, clear description of the algotihm.

Pseudo-arclength continuation

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

Pseudo-arclength continuation was independently developed by Edward Riks and Gerald Wempner for finite element applications in the late 1960s, and published in journals in the early 1970s by H.B. Keller. A detailed account of these early developments is provided in the textbook by M. A. Crisfield: Nonlinear Finite Element Analysis of Solids and Structures, Vol 1: Basic Concepts, Wiley, 1991. Crisfield was one of the most active developers of this class of methods, which are by now standard procedures of commercial nonlinear finite element programs.

PseudoArclength.gif

The algorithm is a predictor-corrector method. The prediction step finds the point (in IR^(n+1) ) which is a step along the tangent vector at the current pointer. The corrector is usually Newton's method, or some variant, to solve the nonlinear system

where is the tangent vector at . The Jacobian of this system is the bordered matrix

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 . 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 became 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). [2]

It is extremely insightful as to the presence of stable solutions (attracting or repelling) in the study of nonlinear differential equations where time 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 Reynolds numbers. Also, 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

Related Research Articles

Quadratic programming (QP) is the process of solving certain mathematical optimization problems involving quadratic functions. Specifically, one seeks to optimize a multivariate quadratic function subject to linear constraints on the variables. Quadratic programming is a type of nonlinear programming.

<span class="mw-page-title-main">Singular value decomposition</span> Matrix decomposition

In linear algebra, the singular value decomposition (SVD) is a factorization of a real or complex matrix into a rotation, followed by a rescaling followed by another rotation. It generalizes the eigendecomposition of a square normal matrix with an orthonormal eigenbasis to any matrix. It is related to the polar decomposition.

<span class="mw-page-title-main">Lyapunov exponent</span> The rate of separation of infinitesimally close trajectories

In mathematics, the Lyapunov exponent or Lyapunov characteristic exponent of a dynamical system is a quantity that characterizes the rate of separation of infinitesimally close trajectories. Quantitatively, two trajectories in phase space with initial separation vector diverge at a rate given by

<span class="mw-page-title-main">Optimal control</span> Mathematical way of attaining a desired output from a dynamic system

Optimal control theory is a branch of control theory that deals with finding a control for a dynamical system over a period of time such that an objective function is optimized. It has numerous applications in science, engineering and operations research. For example, the dynamical system might be a spacecraft with controls corresponding to rocket thrusters, and the objective might be to reach the Moon with minimum fuel expenditure. Or the dynamical system could be a nation's economy, with the objective to minimize unemployment; the controls in this case could be fiscal and monetary policy. A dynamical system may also be introduced to embed operations research problems within the framework of optimal control theory.

In mathematics and computing, the Levenberg–Marquardt algorithm, also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. These minimization problems arise especially in least squares curve fitting. The LMA interpolates between the Gauss–Newton algorithm (GNA) and the method of gradient descent. The LMA is more robust than the GNA, which means that in many cases it finds a solution even if it starts very far off the final minimum. For well-behaved functions and reasonable starting parameters, the LMA tends to be slower than the GNA. LMA can also be viewed as Gauss–Newton using a trust region approach.

Ridge regression is a method of estimating the coefficients of multiple-regression models in scenarios where the independent variables are highly correlated. It has been used in many fields including econometrics, chemistry, and engineering. Also known as Tikhonov regularization, named for Andrey Tikhonov, it is a method of regularization of ill-posed problems. It is particularly useful to mitigate the problem of multicollinearity in linear regression, which commonly occurs in models with large numbers of parameters. In general, the method provides improved efficiency in parameter estimation problems in exchange for a tolerable amount of bias.

In mathematics, linearization is finding the linear approximation to a function at a given point. The linear approximation of a function is the first order Taylor expansion around the point of interest. In the study of dynamical systems, linearization is a method for assessing the local stability of an equilibrium point of a system of nonlinear differential equations or discrete dynamical systems. This method is used in fields such as engineering, physics, economics, and ecology.

<span class="mw-page-title-main">Gauss–Newton algorithm</span> Mathematical algorithm

The Gauss–Newton algorithm is used to solve non-linear least squares problems, which is equivalent to minimizing a sum of squared function values. It is an extension of Newton's method for finding a minimum of a non-linear function. Since a sum of squares must be nonnegative, the algorithm can be viewed as using Newton's method to iteratively approximate zeroes of the components of the sum, and thus minimizing the sum. In this sense, the algorithm is also an effective method for solving overdetermined systems of equations. It has the advantage that second derivatives, which can be challenging to compute, are not required.

In mathematics, a differential-algebraic system of equations (DAE) is a system of equations that either contains differential equations and algebraic equations, or is equivalent to such a system.

In numerical analysis, the Crank–Nicolson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations. It is a second-order method in time. It is implicit in time, can be written as an implicit Runge–Kutta method, and it is numerically stable. The method was developed by John Crank and Phyllis Nicolson in the mid 20th century.

In applied mathematics, in particular the context of nonlinear system analysis, a phase plane is a visual display of certain characteristics of certain kinds of differential equations; a coordinate plane with axes being the values of the two state variables, say (x, y), or (q, p) etc. (any pair of variables). It is a two-dimensional case of the general n-dimensional phase space.

In linear algebra, it is often important to know which vectors have their directions unchanged by a given linear transformation. An eigenvector or characteristic vector is such a vector. Thus an eigenvector of a linear transformation is scaled by a constant factor when the linear transformation is applied to it: . The corresponding eigenvalue, characteristic value, or characteristic root is the multiplying factor .

In the mathematics of evolving systems, the concept of a center manifold was originally developed to determine stability of degenerate equilibria. Subsequently, the concept of center manifolds was realised to be fundamental to mathematical modelling.

In numerical linear algebra, the method of successive over-relaxation (SOR) is a variant of the Gauss–Seidel method for solving a linear system of equations, resulting in faster convergence. A similar method can be used for any slowly converging iterative process.

In computational chemistry, a constraint algorithm is a method for satisfying the Newtonian motion of a rigid body which consists of mass points. A restraint algorithm is used to ensure that the distance between mass points is maintained. The general steps involved are: (i) choose novel unconstrained coordinates, (ii) introduce explicit constraint forces, (iii) minimize constraint forces implicitly by the technique of Lagrange multipliers or projection methods.

Simplicial continuation, or piecewise linear continuation, is a one-parameter continuation method which is well suited to small to medium embedding spaces. The algorithm has been generalized to compute higher-dimensional manifolds by and.

In linear algebra, eigendecomposition is the factorization of a matrix into a canonical form, whereby the matrix is represented in terms of its eigenvalues and eigenvectors. Only diagonalizable matrices can be factorized in this way. When the matrix being factorized is a normal or real symmetric matrix, the decomposition is called "spectral decomposition", derived from the spectral theorem.

In numerical linear algebra, the alternating-direction implicit (ADI) method is an iterative method used to solve Sylvester matrix equations. It is a popular method for solving the large matrix equations that arise in systems theory and control, and can be formulated to construct solutions in a memory-efficient, factored form. It is also used to numerically solve parabolic and elliptic partial differential equations, and is a classic method used for modeling heat conduction and solving the diffusion equation in two or more dimensions. It is an example of an operator splitting method.

Non-linear least squares is the form of least squares analysis used to fit a set of m observations with a model that is non-linear in n unknown parameters (m ≥ n). It is used in some forms of nonlinear regression. The basis of the method is to approximate the model by a linear one and to refine the parameters by successive iterations. There are many similarities to linear least squares, but also some significant differences. In economic theory, the non-linear least squares method is applied in (i) the probit regression, (ii) threshold regression, (iii) smooth regression, (iv) logistic link regression, (v) Box–Cox transformed regressors ().

<span class="mw-page-title-main">Ordinary differential equation</span> Differential equation containing derivatives with respect to only one variable

In mathematics, an ordinary differential equation (ODE) is a differential equation (DE) dependent on only a single independent variable. As with other DE, its unknown(s) consists of one function(s) and involves the derivatives of those functions. The term "ordinary" is used in contrast with partial differential equations (PDEs) which may be with respect to more than one independent variable, and, less commonly, in contrast with stochastic differential equations (SDEs) where the progression is random.

References

  1. Introduction to Numerical Continuation Methods by Eugene L. Allgower and Kurt Georg Colorado State University 1990
  2. Engelnkemper, S.; Gurevich, S. V.; Uecker, H.; Wetzel, D.; Thiele, U. (7 July 2018). Computational modelling of bifurcations and instabilities in fluid dynamics. Springer. pp. 459–501. arXiv: 1808.02321 . doi:10.1007/978-3-319-91494-7_13. ISBN   9783319914930.

Books

[B1] " Introduction to Numerical Continuation Methods ", Eugene L. Allgower and Kurt Georg, SIAM Classics in Applied Mathematics 45. 2003.

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

[B3] "Lyapunov-Schmidt Methods in Nonlinear Analysis and Applications", Nikolay Sidorov, Boris Loginov, Aleksandr Sinitsyn, and Michail Falaleev, Kluwer Academic Publishers, 2002.

[B4] "Methods of Bifurcation Theory", Shui-Nee Chow and Jack K. Hale, Springer-Verlag 1982.

[B5] "Elements of Applied Bifurcation Theory", Yuri A. Kunetsov, Springer-Verlag Applied Mathematical Sciences 112, 1995.

[B6] "Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields", John Guckenheimer and Philip Holmes, Springer-Verlag Applied Mathematical Sciences 42, 1983.

[B7] "Elementary Stability and Bifurcation Theory", Gerard Iooss and Daniel D. Joseph, Springer-Verlag Undergraduate Texts in Mathematics, 1980.

[B8] "Singularity Theory and an Introduction to Catastrophe Theory", Yung-Chen Lu, Springer-Verlag, 1976.

[B9] "Global Bifurcations and Chaos, Analytic Methods", S. Wiggins, Springer-Verlag Applied Mathematical Sciences 73, 1988.

[B10] "Singularities and Groups in Bifurcation Theory, volume I", Martin Golubitsky and David G. Schaeffer, Springer-Verlag Applied Mathematical Sciences 51, 1985.

[B11] "Singularities and Groups in Bifurcation Theory, volume II", Martin Golubitsky, Ian Stewart and David G. Schaeffer, Springer-Verlag Applied Mathematical Sciences 69, 1988.

[B12] "Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems", Alexander Morgan, Prentice-Hall, Englewood Cliffs, N.J. 1987.

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

[B14] "The Implicit Function Theorem: History, Theory and Applications", Steven G. Krantz and Harold R. Parks, Birkhauser, 2002.

[B15] "Nonlinear Functional Analysis", J. T. Schwartz, Gordon and Breach Science Publishers, Notes on Mathematics and its Applications, 1969.

[B16] "Topics in Nonlinear Functional Analysis", Louis Nirenberg (notes by Ralph A. Artino), AMS Courant Lecture Notes in Mathematics 6, 1974.

[B17] "Newton Methods for Nonlinear Problems -- Affine Invariance and Adaptive Algorithms", P. Deuflhard, Series Computational Mathematics 35, Springer, 2006.

Journal articles

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

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

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

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

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

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

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

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

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

[A10] "On the Computation of Multi-Dimensional Solution Manifolds of Parameterized Equations", W.C. Rheinboldt, Numerishe Mathematik, 53, 1988, pages 165-181.

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

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

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

[A14] "Multiple Parameter Continuation: Computing Implicitly Defined k-manifolds", M. E. Henderson, IJBC 12[3]:451-76, 2003.

[A15] "MANPACK: a set of algorithms for computations on implicitly defined manifolds", W. C. Rheinboldt, Comput. Math. Applic. 27 pages 15–9, 1996.

[A16] "CANDYS/QA - A Software System For Qualitative Analysis Of Nonlinear Dynamical Systems", Feudel, U. and W. Jansen, Int. J. Bifurcation and Chaos, vol. 2 no. 4, pp. 773–794, World Scientific, 1992.