Pseudo arclength continuation

This is one of the continuation methods implemented in the package. It is set by the option PALC(tangent = Bordered()) or PALC(tangent = Secant()) in continuation. See also PALC for more information.

For solving

\[\mathbb R^n\ni F(x,p) = 0 \quad\tag{E}\]

using a Newton algorithm, we miss an equation. The simplest way is to select an hyperplane in the space $\mathbb R^n\times \mathbb R$ passing through $(x_0,p_0)$:

\[N(x, p) := \frac{\theta}{n} \langle x - x_0, dx_0\rangle + (1 - \theta)\cdot(p - p_0)\cdot dp_0 - ds = 0\tag{N}\]

with $\theta\in[0,1]$ and where $ds$ is the pseudo arclength (see [Keller]).

Parameter `θ`

The parameter θ in ContinuationPar is very important. It should be tuned for continuation to work properly especially in the case of large problems where the $\langle x - x_0, dx_0\rangle$ component in the constraint might be favored too much. Also, large θs favour p as the corresponding term in the constraint $N$ involves the term $1-θ$.

Predictor

The possible predictors are listed in Predictors - Correctors.

Corrector

The corrector is the newton algorithm for finding the roots $(x,p)$ of

\[\begin{bmatrix} F(x,p) \\ N(x,p)\end{bmatrix} = 0\tag{PALC}\]

Linear Algebra

Norm

First, the option normC continuation specifies the norm used to evaluate the residual in the following way:

\[max(normC(F(x,p)), |N(x,p)|)<tol.\]

It is thus used as a stopping criterion for the corrector. The dot product (resp. norm) used in $N$ and in the (iterative) linear solvers is LinearAlgebra.dot (resp. LinearAlgebra.norm). It can be changed by importing these functions and redefining it. Note that by default, the $\mathcal L^2$ norm is used.

These details are important because the constraint $N$ incorporates the factor length. For some custom type implementing a Vector space, the dot product could already incorporates the length factor in which case you should either redefine the dot product or change $\theta$.

Dot product

In the constraint $N$ above, the scalar product is in fact saved in BifurcationKit.jl as dotp(x,y) -> dot(x,y)/length(y). This is used in the bordered linear solvers associated to PALC. If you want to use your own dot product, you can pass

dotPALC = BK.DotTheta(mydot),

to continuation. Additionally, you may want to provide the linear operator P such that mydot(x,y) = dot(x, A*y), especially if you intend too use the linear solver MatrixBLS. We refer to BifurcationKit.DotTheta for more details.

Linear problem

Pseudo arclength continuation is based on a newton solver applied to the enlarged problem (PALC). We thus need to solve a linear system of size $n+1$ whereas the user passed a linear solver (in ContinuationPar().newton_options) for a system of size $n$.

The linear solver for the linear problem associated to (PALC) is set by the option linear_algo in continuation: it is one of Bordered linear solvers (BLS).

Step size control

Each time the corrector fails, the step size $ds$ is halved. This has the disadvantage of having lost Newton iterations (which costs time) and imposing small steps (which can be slow as well). To prevent this, the step size is controlled internally with the idea of having a constant number of Newton iterations per point. This is in part controlled by the aggressiveness factor a in ContinuationPar.

References

  • Keller

    Keller, Herbert B. Lectures on Numerical Methods in Bifurcation Problems. Springer, 1988