Preservation of Bifurcations of Hamiltonian Boundary Value Problems Under Discretisation

We show that symplectic integrators preserve bifurcations of Hamiltonian boundary value problems and that non-symplectic integrators do not. We provide a universal description of the breaking of umbilic bifurcations by nonysmplectic integrators. We discover extra structure induced from certain types of boundary value problems, including classical Dirichlet problems, that is useful to locate bifurcations. Geodesics connecting two points are an example of a Hamiltonian boundary value problem, and we introduce the jet-RATTLE method, a symplectic integrator that easily computes geodesics and their bifurcations. Finally, we study the periodic pitchfork bifurcation, a codimension-1 bifurcation arising in integrable Hamiltonian systems. It is not preserved by either symplectic or non-symplectic integrators, but in some circumstances symplecticity greatly reduces the error.

1. Truncated miniversal deformations in Arnold’s nomenclature.

2. One can consider different equivalence relations. Our discussion applies to right-left equivalence. It can also be applied to left equivalence (finer) or any coarser notion like contact equivalence.

3. More generally, a generic singularity \(\nabla g(x^1,\ldots ,x^n)\) of the exact flow is broken in the numerical flow of a non-symplectic integrator if and only if for a versal unfolding \((g_\mu )_\mu \) of g the family \((\nabla g_\mu (x^1,\ldots ,x^n))_\mu \) does not constitute a versal roots-of-a-function-type unfolding.

4. Vector fields along a geodesic arising as variational vector fields for variations through geodesics.

5. A 2n-dimensional Hamiltonian system is completely integrable if it possesses n functionally independent, Poisson commuting integrals of motion.

We thank Peter Donelan, Bernd Krauskopf, and Hinke Osinga for useful discussions. This research was supported by the Marsden Fund of the Royal Society Te Apārangi.

Authors and Affiliations

1. School of Fundamental Sciences, Massey University, Palmerston North, New Zealand

   Robert I McLachlan & Christian Offen

Corresponding author

Correspondence to Christian Offen.

Communicated by Hans Munthe-Kaas.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Jet-RATTLE for Calculation of Geodesics on Hypersurfaces

Let (N, g) be a hypersurface of \({\mathbb {R}}^n\) defined by the equation \(f(q)=0\) for \(f:{\mathbb {R}}^n \rightarrow {\mathbb {R}}\) such that \(\nabla f(q) \not =0\) for all \(q \in M\). Here g refers to the induced Riemannian metric on the hypersurface N. In order to compute geodesics on (N, g) with respect to the Levi-Civita connection, we apply the RATTLE method to the Hamiltonian

defined on the cotangent bundle over \(T^*{\mathbb {R}}^n\) with the standard symplectic structure for cotangent bundles and Darboux coordinates \(q^1,\ldots ,q^n,p_1,\ldots ,p_n\). In the above formula, \(\langle .,.\rangle \) denotes the Euclidean scalar product. For a fixed time-step \(h>0\), the RATTLE method gives rise to a map on \(T^*N\) which is symplectic with respect to the standard symplectic structure on cotangent bundles (assuming convergence of the implicit scheme) [17].

The formulas for the time-h-map \(\varPsi _h\) calculating the two n-dimensional vectors \((q_{n+1},p_{n+1})\) from the initial values \((q_{n},p_{n})\) read:

After the 1-dimensional Eq. (29) is solved for \(\lambda \in {\mathbb {R}}\), the remaining equations can be evaluated explicitly.

Remark 9

The formulas (32) and (33) describe a projection of \(p_{n+\frac{1}{2}}\) to the tangent space at \(q_{n+1}\). The effect is wiped out by (29, 30, 31) of the following step, i.e. the value for \(q_{n+2}\) does not depend on whether we set \(p_{n+1}\) according to (33) or simply \(p_{n+1} = p_{n+\frac{1}{2}}\). If the formulas are iterated, then the projection step (32, 33) is only needed in the last step of the iteration (unless one is interested in the intermediate values for p themselves). Indeed, in the examples presented in this paper not only the intermediate p-values but also the final momentum is irrelevant. This means for the calculation of conjugate loci one could simply use

The derivative \(D\varPsi _h\) of the time-h-map (including the projection step) can be obtained by evaluating the following formulas. We interpret the vectors \(q_n\), \(p_n\) and the gradient vectors \(\nabla f(q_n)\), \(\nabla _q \lambda \), etc. as column vectors such that, for instance, \(\nabla f(q_n) (\nabla _q \lambda )^T\) denotes a dyadic product. The symbol I refers to an n-dimensional identity matrix.

When the time-h-map \(\varPsi _h\) is iterated N-times to obtain the numerical time-Nh-map \(\varPhi \), the derivatives can be updated as follows:

We obtain the derivatives as \(D\varPhi (q_0,p_0)=V^{(N)}\). We refer to this 1-jet version of the RATTLE method applied to a hypersurface as jet-RATTLE.

Breaking of an Elliptic Umbilic Using a Non-symplectic Integrator: Numerical Example

Let us compare the capturing of an elliptic umbilic bifurcation \(D_4^-\) by the second-order accurate symplectic Störmer–Verlet method to a non-symplectic method of the same order of accuracy. For this, we consider the Dirichlet problem for the Hénon–Heiles-type Hamiltonian system described in Sect. 3.3. In contrast to the numerical experiment described in 3.3, we reduce the number of time-steps from \(N=10\) to \(N=5\) and perturb the Hamiltonian from (20) with the extra term \(0.01 y_2 \sin (y_1)\) to

Resolving an elliptic umbilic bifurcation \(D^-_4\) with the symplectic Störmer–Verlet method. The plot to the left shows the set B, and the plot to the right shows the level bifurcation set \({\mathcal {B}}\)

Full size image

Resolving an elliptic umbilic bifurcation \(D^-_4\) with the non-symplectic second-order Runge–Kutta method. The plot to the left shows the set B and the plot to the right shows the level bifurcation set \({\mathcal {B}}\). The set was rotated around the \(Y_2\) axis by \(0.0271 \mathrm {rad}\) in order to allow for a convenient rescaling of the axes. Instead of an elliptic umbilic bifurcation, there are three lines of cusp bifurcations which fail to merge

Full size image

In the considered boundary value problem

where \(\phi ^X=(\phi ^{X^1},\phi ^{X^2})\) are the x-components of the Hamiltonian flow map at time 1 and \((x^*)^2\), \((X^*)^1\), \((X^*)^2\) the parameters of the problem; the level bifurcation set is locally given by

for a subset \(U \subset {\mathbb {R}}^4\) of the phase space. The level bifurcation set \({\mathcal {B}}\) can be obtained from

using \(\phi ^X\). Figures 23 and 24 show plots of the sets B (to the left) and \({\mathcal {B}}\) (to the right). For Fig. 23, the flow \(\phi \) was approximated with the symplectic Störmer–Verlet method. We see an elliptic umbilic bifurcation, where three lines of cusps merge in one singular point marked by an asterisk. Its position in the phase portrait of the numerical flow can be calculated as a root of \((x^2,y_1,y_2) \mapsto D_y \phi ^X(0,x^2,y_1,y_2)\). For Fig. 24, the flow \(\phi \) was approximated with the explicit midpoint rule (RK2), which is a second-order non-symplectic Runge–Kutta method. While the sheets in the plot for B still approach a singular point, they cannot reach it and connect in a whole circle rather than a singular point. In the level bifurcation set, this has the effect that we do not obtain an elliptic umbilic bifurcation but three lines of cusp bifurcations which fail to merge in an umbilic point.

To know which parts of the bifurcation diagrams are to be compared, the calculations were first done to high accuracy such that the bifurcation diagrams obtained by the Störmer–Verlet method and by RK2 were close. Then, the step-sizes were increased gradually and the movement of the singular point where the matrix \(D_y \phi ^X(0,x^2,y_1,y_2)\) is near the zero matrix was tracked in both simulations.

Reprints and permissions