Numerical continuation and bifurcation

Let an algebraic problem coming from discretisation of an FEM-model can be written in the form

\[F(U) = 0.\]

In what follows, we shall suppose that the model depends on an additional scalar parameter \(\lambda\) so that \(F(U) = F(U, \lambda)\).

Numerical continuation

Methods of numerical continuation serve for tracing solutions of the system

\[F(U, \lambda) = 0, \quad F\colon \mathbb{R}^{N} \times \mathbb{R} \to \mathbb{R}^{N}.\]

In GetFEM, a continuation technique for piecewise \(C^{1}\) (\(PC^{1}\)) solution curves is implemented (see [Li-Re2014] for more details). Since it does not make an explicit difference between the state variable \(U\) and the parameter \(\lambda\), we shall denote \(Y := (U, \lambda)\) for brevity. Nevertheless, to avoid bad scaling when calculating tangents, for example, we shall use the following weighted scalar product and norm:

\[\langle Y, \tilde{Y} \rangle_{w} := \kappa \langle U, \tilde{U} \rangle + \lambda \tilde{\lambda},\quad \lVert Y \rVert_{w} := \sqrt{\kappa \lVert U \rVert^{2} + \lambda^{2}},\qquad Y = (U, \lambda),\, \tilde{Y} = (\tilde{U}, \tilde{\lambda}).\]

Here, \(\kappa\) should be chosen so that \(\kappa \langle U, \tilde{U} \rangle\) is proportional to the scalar product of the corresponding space variables, usually in \(L^{2}\). One can take, for example, \(\kappa = h^{d}\), where \(h\) is the mesh size and \(d\) stands for the dimension of the underlying problem. Alternatively, \(\kappa\) can be chosen as \(1/N\) for simplicity.

The idea of the continuation strategy is to continue smooth pieces of solution curves by a classical predictor-corrector method and to join the smooth pieces continuously.

The particular predictor-corrector method employed is a slight modification of the inexact Moore-Penrose continuation implemented in MATCONT [Dh-Go-Ku2003]. It computes a sequence of consecutive points \(Y_{j}\) lying approximately on a solution curve and a sequence of the corresponding unit tangent vectors \(T_{j}\):

\[\lVert F(Y_{j}) \rVert \leq \varepsilon,\quad F'(Y_{j}; T_{j}) = 0,\quad \lVert T_{j} \rVert_{w} = 1,\quad j = 0, 1,\dotsc.\]

To describe it, let us suppose that we have a couple \((Y_{j}, T_{j})\) satisfying the relations above at our disposal. In the prediction, an initial approximation of \((Y_{j+1}, T_{j+1})\) is taken as

\[Y_{j+1}^{0} := Y_{j} + h_{j} T_{j},\quad T_{j+1}^{0} := T_{j},\]

where \(h_{j}\) is a step size. Its choice will be discussed later on.

In the correction, one computes a sequence \(\{(Y_{j+1}^{l}, T_{j+1}^{l})\}\), where \(T_{j+1}^{l} := \tilde{T}_{j+1}^{l} / \lVert \tilde{T}_{j+1}^{l} \rVert_{w}\) and the couple \((Y_{j+1}^{l}, \tilde{T}_{j+1}^{l})\) is given by one iteration of the Newton method applied to the equation \(F^{l}(Y, T) = 0\) with

\[\begin{split}F^{l}(Y, T) := \begin{pmatrix}F(Y)\\ (T_{j+1}^{l-1})^{\top}(Y - Y_{j+1}^{l-1})\\ \nabla F(Y_{j+1}^{l-1})T\\ \langle T_{j+1}^{l-1}, T \rangle_{w} - \langle T_{j+1}^{l-1}, T_{j+1}^{l-1} \rangle_{w}\end{pmatrix}\end{split}\]

and the initial approximation \((Y_{j+1}^{l-1}, T_{j+1}^{l-1})\). Due to the potential non-differentiability of \(F\), a piecewise-smooth variant of the Newton method is used (Algorithm 7.2.14 in [Fa-Pa2003]).

Correction.

A couple \((Y_{j+1}^{l}, T_{j+1}^{l})\) is accepted for \((Y_{j+1}, T_{j+1})\) if \(\lVert F(Y_{j+1}^{l})\rVert \leq \varepsilon\), \(\lVert Y_{j+1}^{l} - Y_{j+1}^{l-1}\rVert_{w} \leq \varepsilon'\), and the cosine of the angle between \(T_{j+1}^{l}\) and \(T_{j}\) is greater or equal to \(c_{\mathrm{min}}\). Let us note that the partial gradient of \(F\) (or of one of its selection functions in the case of the non-differentiability) with respect to \(U\) is assembled analytically whereas the partial gradient with respect to \(\lambda\) is evaluated by forward finite differences with an increment equal to 1e-8.

The step size \(h_{j+1}\) in the next prediction depends on how the Newton correction has been successful. Denoting the number of iterations needed by \(l_{\mathrm{it}}\), it is selected as

\[\begin{split}h_{j+1} := \begin{cases}\max\{h_{\mathrm{dec}} h_{j}, h_{\mathrm{min}}\}& \text{if no new couple has been accepted},\\ \min\{h_{\mathrm{inc}} h_{j}, h_{\mathrm{max}}\}& \text{if a new couple has been accepted and } l_{\mathrm{it}} < l_{\mathrm{thr}},\\ h_{j}& \text{otherwise},\end{cases}\end{split}\]

where \(0 < h_{\mathrm{dec}} < 1 < h_{\mathrm{inc}}\), \(0 < l_{\mathrm{thr}}\) and \(0 < h_{\mathrm{min}} < h_{\mathrm{max}}\) are given constants. At the beginning, one sets \(h_{1} := h_{\mathrm{init}}\) for some \(h_{\mathrm{min}} \leq h_{\mathrm{init}} \leq h_{\mathrm{max}}\).

Now, let us suppose that we have approximated a piece of a solution curve corresponding to one sub-domain of smooth behaviour of \(F\) and we want to recover a piece corresponding to another sub-domain of smooth behaviour. Let \((Y_{j},T_{j})\) be the last computed couple.

Transition between smooth pieces of a solution curve.

To approximate the tangent to the other smooth piece, we first take a point \(Y_{j} + h T_{j}\) with \(h\) a bit greater than \(h_{\mathrm{min}}\) so that this point belongs to the interior of the other sub-domain of smooth behaviour. Then we find \(\tilde{T}\) such that

\[\nabla F(Y_{j} + h T_{j}) \tilde{T} = 0,\quad \lVert \tilde{T} \rVert_{w} = 1,\]

and it remains to determine an appropriate direction of this vector. This can be done on the basis of the following observations: First, there exists \(r \in \{\pm 1\}\) such that \(Y_{j} - r \tilde{h} \tilde{T}\) remains in the same sub-domain as \(Y_{j}\) for any \(\tilde{h}\) positive. This is characterised by the fact that \(\frac{\lvert T_{-}^{\top} \tilde{T}\rvert}{\lVert T_{-} \rVert \lVert \tilde{T} \rVert}\) is significantly smaller than 1 for \(T_{-}\) with \(\nabla F(Y_{j} - r \tilde{h} \tilde{T}) T_{-} = 0\). Second, \(Y_{j} + r \tilde{h} \tilde{T}\) appears in the other sub-domain for \(\tilde{h}\) larger than some positive threshold, and, for such values, \(\frac{\lvert T_{+}^{\top} \tilde{T}\rvert}{\lVert T_{+} \rVert \lVert \tilde{T} \rVert}\) is close to 1 for \(T_{+}\) with \(\nabla F(Y_{j} + r \tilde{h} \tilde{T}) T_{+} = 0\).

This suggests the following procedure for selecting the desired direction of \(\tilde{T}\): Increase the values of \(\tilde{h}\) successively from \(h_{\mathrm{min}}\), and when you arrive at \(\tilde{h}\) and \(r \in \{\pm 1\}\) such that

\[\frac{\lvert T^{\top} \tilde{T}\rvert}{\lVert T \rVert \lVert \tilde{T} \rVert} \approx 1\quad \text{if}\ \nabla F(Y_{j} + r \tilde{h} \tilde{T}) T = 0,\]

take \(r \tilde{T}\) as the approximation of the tangent to the other smooth piece.

Having this approximation at our disposal, we restart the predictor-corrector with \((Y_{j}, r \tilde{T})\).

In GetFEM, the continuation is implemented for two ways of parameterization of the model:

1. The parameter \(\lambda\) is directly a scalar datum, which the model depends on.

2. The model is parametrised by the scalar parameter \(\lambda\) via a vector datum \(P\), which the model depends on. In this case, one takes the linear path

   \[\lambda \mapsto P(\lambda) := (1 - \lambda)P^{0} + \lambda P^{1},\]

   where \(P^{0}\) and \(P^{1}\) are given values of \(P\), and one traces the solution set of the problem

   \[F(U, P(\lambda)) = 0.\]

Detection of limit points

When tracing solutions of the system \(F(U,\lambda) = 0\), one may be interested in limit points (also called fold or turning points), where the number of solutions with the same value of \(\lambda\) changes. These points can be detected by a sign change of a test function \(\tau_{\mathrm{LP}}\):

\[\tau_{\mathrm{LP}}(T_{j}) \tau_{\mathrm{LP}}(T_{j+1}) < 0,\]

where \(\tau_{\mathrm{LP}}\) is defined by

\[\tau_{\mathrm{LP}}(T) := T_{\lambda},\quad T = (T_{U},T_{\lambda}) \in \mathbb{R}^{N} \times \mathbb{R}.\]

Limit point.

Numerical bifurcation

A point \(\bar{Y}\) is called a bifurcation point of the system \(F(Y) = 0\) if \(F(\bar{Y}) = 0\) and two or more distinct solution curves pass through it. The following result gives a test for smooth bifurcation points (see, e.g., [Georg2001]):

Let \(s \mapsto Y(s)\) be a parameterization of a solution curve and \(\bar{Y} := Y(\bar{s})\) be a bifurcation point. Moreover, let \(T^{\top} \dot{Y}(\bar{s}) > 0\), \(B \notin \mathrm{Im}(J(\bar{Y}))\), \(C \notin \mathrm{Im}(J(\bar{Y})^{\top})\), \(d \in \mathbb{R}\) and

\[\begin{split}J(Y) := \begin{pmatrix}\nabla F(Y)\\ T^{\top}\end{pmatrix}.\end{split}\]

Define \(\tau_{\mathrm{BP}}(Y)\) via

\[\begin{split}\begin{pmatrix}J(Y)& B\\ C^{\top}& d\end{pmatrix} \begin{pmatrix}V(Y)\\ \tau_{\mathrm{BP}}(Y)\end{pmatrix} = \begin{pmatrix}0\\ 1\end{pmatrix}.\end{split}\]

Then \(\tau_{\mathrm{BP}}(Y(s))\) changes its sign at \(s = \bar{s}\).

Obviously, if one takes \(B\), \(C\) and \(d\) randomly, it is highly possible that they satisfy the requirements above. Consequently, the numerical continuation method is able to detect bifurcation points by taking the vectors \(Y\) and \(T\) supplied by the correction at each continuation step and monitoring the signs of \(\tau_{\mathrm{BP}}\).

Once a bifurcation point \(\bar{Y}\) is detected by a sign change \(\tau_{\mathrm{BP}}(Y_{j}) \tau_{\mathrm{BP}}(Y_{j+1}) < 0\), it can be approximated more precisely by the predictor-corrector steps described above with a special step-length adaptation (see Section 8.1 in [Al-Ge1997]). Namely, one can take the subsequent step lengths as

\[h_{j+1} := -\frac{\tau_{\mathrm{BP}}(Y_{j+1})}{\tau_{\mathrm{BP}}(Y_{j+1}) - \tau_{\mathrm{BP}}(Y_{j})}h_{j}\]

until \(\lvert h_{j+1} \rvert < h_{\mathrm{min}}\), which corresponds to the secant method for finding a zero of the function \(s \mapsto \tau_{\mathrm{BP}}(Y(s))\).

Finally, it would be desirable to switch solution branches. To this end, we shall consider the case of the so-called simple bifurcation point, where only two distinct solution curves intersect.

Let \(\tilde{Y}\) be an approximation of \(\bar{Y}\) that we are given and \(V(\tilde{Y})\) be the first part of the solution of the augmented system for computing the test function \(\tau_{\mathrm{BP}}(\tilde{Y})\). As proposed in [Georg2001], one can take \(V(\tilde{Y})\) as a predictor direction and do one continuation step starting with \((\tilde{Y}, V(\tilde{Y}))\) to obtain a point on a new branch. After this continuation step has been performed successfully and a point on the new branch has been recovered, one can proceed with usual predictor-corrector steps to trace this branch.

Recently, tools for numerical \(PC^{1}\)-bifurcation have been developed in GetFEM. Let \(J\) be a matrix function of a real parameter now defined by

\[\begin{split}J(\alpha) := (1-\alpha)\begin{pmatrix}\nabla F(Y_{j})\\ T_{j}^{\top}\end{pmatrix} + \alpha\begin{pmatrix}\nabla F(Y_{j+1})\\ T_{j+1}^{\top}\end{pmatrix}.\end{split}\]

As proposed in [Li-Re2014hal], the following test can be used for detection of a \(PC^{1}\) bifurcation point between \(Y_{j}\) and \(Y_{j+1}\):

\[\det J(0) \det J(1) < 0.\]

To perform this test numerically, introduce

\[\begin{split}M(\alpha) := \begin{pmatrix}J(\alpha)& B\\ C^{\top}& d\end{pmatrix}\end{split}\]

and \(\tau_{\mathrm{BP}}(\alpha)\) analogously as above via

\[\begin{split}M(\alpha) \begin{pmatrix}V(\alpha)\\ \tau_{\mathrm{BP}}(\alpha)\end{pmatrix} = \begin{pmatrix}0\\ 1\end{pmatrix}.\end{split}\]

It follows from Cramer’s rule that

\[\tau_{\mathrm{BP}}(\alpha) = \frac{\det J(\alpha)}{\det M(\alpha)}\]

provided that \(\det M(\alpha)\) is non-zero. Hence if \(B\), \(C\) and \(d\) are chosen so that \(\det M(\alpha)\) is non-zero whenever \(\det J(\alpha)\) is zero, then the sign changes of \(\det J(\alpha)\) are characterised by passings of \(\tau_{\mathrm{BP}}(\alpha)\) through 0 whereas the sign changes of \(\det M(\alpha)\) by sign changes of \(\tau_{\mathrm{BP}}(\alpha)\) caused by singularities. To conclude, the sign of \(\det J(0)\det J(1)\) is determined by following the behaviour of \(\tau_{\mathrm{BP}}(\alpha)\) and monitoring the sign changes of \(\det J(\alpha)\) when \(\alpha\) passes through \([0,1]\).

As justified in [Li-Re2014hal], \(B\), \(C\) and \(d\) can be chosen randomly again. The increments \(\delta\) of the current values of \(\alpha\) are changed adaptively so that singularities of \(\tau_{\mathrm{BP}}\) are treated effectively. After each calculation of \(\tau_{\mathrm{BP}}(\alpha)\), \(\delta\) is set as follows:

\[\begin{split}\delta := \begin{cases}\min\{2\delta, \delta_{\mathrm{max}}\}& \text{if $\lvert \tau_{\mathrm{BP}}(\alpha) - \tau_{\mathrm{BP}}(\alpha - \delta) \rvert < 0.5 \tau_{\mathrm{fac}} \tau_{\mathrm{ref}}$,}\\ \max\{0.1\delta, \delta_{\mathrm{min}}\}& \text{if $\lvert \tau_{\mathrm{BP}}(\alpha) - \tau_{\mathrm{BP}}(\alpha - \delta) \rvert > \tau_{\mathrm{fac}} \tau_{\mathrm{ref}}$,}\\ \delta& \text{otherwise},\end{cases}\end{split}\]

where \(\delta_{\mathrm{max}} > \delta_{\mathrm{min}} > 0\) and \(\tau_{\mathrm{fac}} > 0\) are given constants and \(\tau_{\mathrm{ref}} := \max\{\lvert \tau_{\mathrm{BP}}(1) - \tau_{\mathrm{BP}}(0) \rvert, 10^{-8}\}\).

When a \(PC^{1}\) bifurcation point is detected between \(Y_{j}\) and \(Y_{j+1}\), it is approximated more precisely by a bisection-like procedure. The obtained approximation lies on the same smooth branch as \(Y_{j},\) and the corresponding unit tangent that points out from the corresponding region of smoothness is calculated too.

Contrary to the smooth case, it is not clear how many branches can emanate from the \(PC^{1}\) bifurcation point and in which directions they could be sought. For this reason, continuation steps for a whole sequence of predictor directions are tried out for finding points on new branches.

Denoting \(\tilde{Y}\), \(\tilde{T}\) the approximation of the bifurcation point and the corresponding tangent, respectively, the predictor directions are taken as follows: For a couple of reference vectors \(\tilde{V}_{1}\) and \(\tilde{V}_{2}\), one takes \(\pm V\) with \(V\) satisfying

\[\nabla F(\tilde{Y}+h_{\mathrm{min}}\tilde{V}) V = 0, \quad \lVert V \rVert_{w} = 1,\]

where \(\tilde{V}\) passes through a set of linear combinations of \(\tilde{V}_{1}\) and \(\tilde{V}_{2}\). The total number of the linear combinations is given by \(n_{\mathrm{dir}},\) and the reference vectors are chosen successively according to the following strategy:

1. One takes \(\tilde{V}_{1} := -\tilde{T}\) and \(\tilde{V}_{2}\) such that

   \[\nabla F(\tilde{Y}+h_{\mathrm{min}}\tilde{T}) \tilde{V}_{2} = 0, \quad \lVert \tilde{V}_{2} \rVert_{w} = 1.\]

2. Let \(\{\tilde{T}_{1},\dotsc\tilde{T}_{n_{\mathrm{br}}}\}\) denote the set of unit tangents that correspond to the points from the branches found so far and that are oriented in the directions of branching from the bifurcation point. Then \(\tilde{V}_{1}\) and \(\tilde{V}_{2}\) are taken successively as different combinations from \(\{\tilde{T}_{1},\dotsc\tilde{T}_{n_{\mathrm{br}}}\}\).

3. If all combinations that are available so far have already been used, let \(\tilde{V}_{1}\) be unchanged and take \(\tilde{V}_{2} := \tilde{V}_{2}^{+}\) with \(\tilde{V}_{2}^{+}\) satisfying

   \[\nabla F\Bigl(\tilde{Y}+h_{\mathrm{min}}\Bigl(\tilde{V}_{2}^{-} + 0.1\frac{\tilde{V}_{3}}{\lVert \tilde{V}_{3} \rVert_{w}}\Bigr)\Bigr) \tilde{V}_{2}^{+} = 0, \quad \lVert \tilde{V}_{2}^{+} \rVert_{w} = 1.\]

   Here, \(\tilde{V}_{2}^{-}\) equals the vector \(\tilde{V}_{2}\) employed previously and \(\tilde{V}_{3}\) is chosen randomly.

The total number of selections of \(\tilde{V}_{1}\) and \(\tilde{V}_{2}\) is given by \(n_{\mathrm{span}}\).

More details on \(PC^1\) numerical branching can be found in [Li-Re2015hal].

Approximation of solution curves of a model

The numerical continuation is defined in getfem/getfem_continuation.h. In order to use it, one has to set it up via the corresponding object first:

getfem::cont_struct_getfem_model S(model, parameter_name, sfac, ls, h_init, h_max, h_min, h_inc, h_dec,
                                   maxit, thrit, maxres, maxdiff, mincos, maxres_solve, noisy, singularities,
                                   non-smooth, delta_max, delta_min, thrvar, ndir, nspan);

where parameter_name is the name of the model datum representing \(\lambda\), sfac represents the scale factor \(\kappa\), and ls is the name of the solver to be used for the linear systems incorporated in the process (e.g., getfem::default_linear_solver<getfem::model_real_sparse_matrix, getfem::model_real_plain_vector>(model)). The real numbers h_init, h_max, h_min, h_inc, h_dec denote \(h_{\mathrm{init}}\), \(h_{\mathrm{max}}\), \(h_{\mathrm{min}}\), \(h_{\mathrm{inc}}\), and \(h_{\mathrm{dec}}\), the integer maxit is the maximum number of iterations allowed in the correction and thrit, maxres, maxdiff, mincos, and maxres_solve denote \(l_{\mathrm{thr}}\), \(\varepsilon\), \(\varepsilon'\), \(c_{\mathrm{min}}\), and the target residual value for the linear systems to be solved, respectively. The non-negative integer noisy determines how detailed information has to be displayed in the course of the continuation process (the larger value the more details), the integer singularities determines whether the tools for detection and treatment of singular points have to be used (0 for ignoring them completely, 1 for detecting limit points, and 2 for detecting and treating bifurcation points, as well), and the boolean value of non-smooth determines whether only tools for smooth continuation and bifurcation have to be used or even tools for non-smooth ones do. The real numbers delta_max, delta_min and thrvar represent \(\delta_{\mathrm{max}}\), \(\delta_{\mathrm{min}}\) and \(\tau_{\mathrm{fac}}\), and the integers ndir and nspan stand for \(n_{\mathrm{dir}}\) and \(n_{\mathrm{span}}\), respectively.

Optionally, parameterization by a vector datum is then declared by:

S.set_parametrised_data_names(initdata_name, finaldata_name, currentdata_name);

Here, the data names initdata_name and finaldata_name should represent \(P^{0}\) and \(P^{1}\), respectively. Under currentdata_name, the values of \(P(\lambda)\) have to be stored, that is, actual values of the datum the model depends on.

Next, the continuation is initialised by:

S.init_Moore_Penrose_continuation(U, lambda, T_U, T_lambda, h);

where U should be a solution for the value of the parameter \(\lambda\) equal to lambda so that \(Y_{0}=\) (U,lambda). During this initialisation, an initial unit tangent \(T_{0}\) corresponding to \(Y_{0}\) is computed in accordance with the sign of the initial value T_lambda, and it is returned in T_U, T_lambda. Moreover, h is set to the initial step size h_init.

Subsequently, one step of the continuation is called by

S.Moore_Penrose_continuation(U, lambda, T_U, T_lambda, h, h0);

After each call, a new point on a solution curve and the corresponding tangent are returned in the variables U, lambda and T_U, T_lambda. The step size for the next prediction is returned in h. The size of the current step is returned in the optional argument h0. According to the chosen value of singularities, the test functions for limit and bifurcation points are evaluated at the end of each continuation step. Furthermore, if a smooth bifurcation point is detected, the procedure for numerical bifurcation is performed and an approximation of the branching point as well as tangents to both bifurcating curves are saved in the continuation object S. From there, they can easily be recovered with member functions of S so that one can initialise the continuation to trace either of the curves next time.

Complete examples of use on a smooth problem are shown in the test programs tests/test_continuation.cc, interface/tests/matlab/demo_continuation.m and interface/src/scilab/demos/demo_continuation.sce, whereas interface/src/scilab/demos/demo_continuation_vee.sce and interface/src/scilab/demos/demo_continuation_block.sce employ also non-smooth tools.