3.11.6. Jumping on bifurcations

If one is just interested at the critical parameter value where the bifurcation happens, e.g. \(r=0\) in all previously discussed cases, there is another possibility to spot the bifurcation than performing arclength continuation and the calculation of eigenvalues. In complicated systems, one neither knows the critical stationary solution \(\vec{x}_\text{c}\) (here a vector comprising e.g. multiple coupled ODEs) nor the critical parameter \(r_\text{c}\). However, one knows that the following has to hold at e.g. a fold bifurcation:

\[\begin{split}\begin{aligned} \vec{F}(\vec{x}_\text{c},r_\text{c})&=0 \quad \text{i.e. }\vec{x}_\text{c}\text{ is a stationary solution at the parameter }r_\text{c}\\ \mathbf{J}(\vec{x}_\text{c},r_\text{c})\vec{v}&=0 \quad \text{i.e. there is an eigenvector }\vec{v}\text{ with eigenvalue zero} \end{aligned}\end{split}\]

For a system with \(N\) degrees of freedom (\(N=\operatorname{dim}(\vec{x})\)), we have now \(2N+1\) unknowns, namely the \(N\) unknowns of \(\vec{x}_\text{c}\), the \(N\) components of the unknown eigenvector \(\vec{v}\) and the critical parameter \(r_\text{c}\), for which an eigenvalue becomes zero. There is obviously one equations missing, which related with the magnitude of \(\vec{v}\), since without any further equation all parameters \(r\) with corresponding stationary solutions will solve it for the trivial choice \(\vec{v}=0\). One could demand that e.g. \(\|\vec{v}\|=1\), but since a reasonable initial guess \(\vec{v}_\text{g}\) for the eigenvector \(\vec{v}\) is usually required anyhow, we just demand \(\vec{v}\cdot\vec{v}_\text{g}=1\). Thereby, this \((2N+1)^\text{th}\) equation is linear.

For the fold normal form (3.14), one gets the system

\[\begin{aligned} r_\text{c}-x^2_\text{c}=0, \qquad -2x_\text{c}v=0, \qquad vv_\text{g}=1 \end{aligned}\]

which yields for all \(v_\text{g}\neq 0\) the known solution \(r_\text{c}=0\) and \(x_\text{c}=0\). For more complicated systems, however, this has to be solved numerically, which can be done in pyoomph by the activate_bifurcation_tracking() method (again after importing the problem and equation classes from bifurcation_fold_param_change.py):

from bifurcation_fold_param_change import *

if __name__=="__main__":
    with FoldProblem() as problem:

        # Find any start solution, which must be close to the bifurcation
        problem.r.value=1
        problem.get_ode("fold").set_value(x=1)
        problem.solve()

        # Find a guess for the normalization constraint
        problem.solve_eigenproblem(0)
        vguess=problem.get_last_eigenvectors()[0] # use the eigenvector as guess

        # Activate fold bifurcation tracking in parameter r and solve the augmented system
        problem.activate_bifurcation_tracking(problem.r,"fold",eigenvector=vguess)
        problem.solve()

        print(f"Critical at r_c={problem.r.value} and x_c={problem.get_ode('fold').get_value('x')}")

Download this example

Download all examples

The same works also for a pitchfork bifurcation, but these are subject to a symmetry, which is broken by the bifurcation. If we apply the fold tracking method to the pitchfork normal form (3.16), we would get the augmented system

\[\begin{aligned} r_\text{c}x_\text{c}\pm x^3_\text{c}=0, \qquad \left(r\pm 3x^2_\text{c}\right)v=0, \qquad vv_\text{g}=1\,. \end{aligned}\]

While it indeed gives the correct solution \(x_\text{c}=0\) at \(r_\text{c}=0\), the root has a multiplicity of three. Numerically, this hampers the Newton solver to converge.

To reduce the multiplicity and account for the symmetry in the pitchfork bifurcation, a symmetry vector \(\vec{\psi}\) is considered and following the \(2N+2\)-system is solved:

\[\begin{split}\begin{aligned} \vec{F}(\vec{x}_\text{c},r_\text{c})+\sigma\vec{\psi}&=0 \quad \text{i.e. }\vec{x}_\text{c}\text{ is a stationary solution at the parameter }r_\text{c}\\ \mathbf{J}(\vec{x}_\text{c},r_\text{c})\vec{v}&=0 \quad \text{i.e. there is an eigenvector }\vec{v}\text{ with eigenvalue zero} \\ \vec{v}_\text{g}\cdot\vec{v}&=1 \quad \text{i.e. the eigenvector }\vec{v}\text{ in non-trivial (for reasonables guesses) }\vec{v}_\text{g}\\ \vec{\psi}\cdot\vec{x}_\text{c}&=0 \quad \text{i.e. the solution is symmetric with respect to the symmetry vector }\vec{\psi} \end{aligned}\end{split}\]

Note that the slack variable \(\sigma\) enforcing the symmetry will be zero at the bifurcation. For the pitchfork normal form with the simple scalar equivalents \(\psi=v_\text{g}=1\), we obtain indeed \(x_\text{c}=0,\,r_\text{c}=0,\,v=1,\,\sigma=0\) and the root of the system has a multiplicity of one, i.e. the Jacobian of the augmented system at the solution is invertible. Thereby, the Newton solver converges well. To find a pitchfork bifurcation, one just has to pass "pitchfork" instead of "fold" as second argument for the activate_bifurcation_tracking() method. Again, one can pass an eigenvector argument which will be used as eigenvector normalization vector \(\vec{v}_\text{g}\) and as symmetry vector \(\psi\). Please refer to the supplied code bifurcation_pitchfork_tracking.py (dependent on bifurcation_pitchfork_arclength_eigen.py) for an example.

Warning

When solving with activate_bifurcation_tracking(), you must deactivate it via deactivate_bifurcation_tracking() after the solve to solve the normal system again (i.e. the system without the augmentation). Also the calculation of eigenvalues and -vectors does not work as usual while bifurcation tracking is active. See next section how to obtain the critical eigenvector.

Warning

Tracking pitchfork bifurcations in spatio-temporal problems (cf. Section 5) requires that the mesh is conforming with the symmetry vector, i.e. the mesh should be also symmetric along the symmetry that it broken by the bifurcation.

Note

Pyoomph can improve the convergence of bifurcation tracking by calling the method setup_for_stability_analysis(). This will generate symbolical C code for the required Hessian terms in the augmented systems for bifurcation tracking. Without calling this method, finite differences are used to calculate the Hessian terms, which can be less accurate and slower. For more details, we refer to Section 10.2.1 and our article [16].