Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm and long time performance

Explicit symplectic integrators have been important tools for accurate and efficient approximations of mechanical systems with separable Hamiltonians. This article proposes for arbitrary Hamiltonians similar integrators, which are explicit, of any even order, symplectic in an extended phase space, and with pleasant long time properties. They are based on a mechanical restraint that binds two copies of phase space together. Using backward error analysis, Kolmogorov-Arnold-Moser theory, and additional multiscale analysis, an error bound of $O\left(T{\delta }^l\omega \right)$ is established for integrable systems, where $T,\delta ,l$, and $\omega$ are, respectively, the (long) simulation time, step size, integrator order, and some binding constant. For nonintegrable systems with positive Lyapunov exponents, such an error bound is generally impossible, but satisfactory statistical behaviors were observed in a numerical experiment with a nonlinear Schrödinger equation.

Figure 1

Comparison between the proposed integrator and a classical nonsymplectic method. Time step $\delta =0.1$ for both methods, $Q\left(0\right)=-3,P\left(0\right)=0,\omega =20$, and $T=1000$.

Figure 2

Comparison between a low-order version of the proposed integrator and the pioneer nonsymplectic method recommended in [34]. Time step $\delta =0.1$ for both methods, $Q\left(0\right)=-3,P\left(0\right)=0,\omega =20$, and $T=10$.

Figure 3

Comparison between the proposed integrators and the pioneer nonsymplectic method: Schwarzschild geodesic and its dissipated version $\left(\gamma ={10}^{-4}\right)$. Time step $\delta =0.2$ for both methods, $\omega =2,T=50000,M=10,m=1$, initial condition $Q\left(0\right)=[0,20,0],P\left(0\right)\approx [0.982,0,-4.472]$ corresponds to initial semimajor axis $a\left(0\right)=20$, eccentricity $e\left(0\right)=0,p_r\left(0\right)=0,p_{\varphi }\left(0\right)=\sqrt{r\left(0\right)}$, and $p_t\left(0\right)$ solves $H\left(0\right)=m^2/2$. Plotted were maxima of relative errors up to given times, respectively scaled by dividing over the Keplerian period $2\pi \sqrt{a{\left(0\right)}^3/M}$ for $t$, over $a\left(0\right)[1+e(0\left)\right]$ for $r$, over $2\pi$ for $\varphi$, and unscaled for $H$.

Figure 4

Long time simulations of NLS with $N=2$ by fourth-order Runge-Kutta and a proposed integrator ${\varphi }_4^{\delta }$. $[q_1,p_1,q_2,p_2]\left(0\right)=[3,1,0.01,0],\delta =0.01,\omega =100$. Long time convergence of the trajectory towards the ergodic limit is quantified in the right column based on $\langle I_1\rangle \left(t\right)-\langle I_2\rangle \left(t\right)$ [see Eq. (5)]; this observable was chosen to reduce the interference with the numerical loss of total mass in RK4. Computation done using win32 matlab R2010a on x64 Windows 7 i7-4600U CPU.

Figure 5

Short time simulation of NLS with $N=5$ by ${\varphi }_4$. $[q_1,p_1]\left(0\right)=[3,1],[q_i,p_i]\left(0\right)=[0.01,0]$ for $i\ne 1,\delta =0.001$ for a closer examination of details, $\omega =100$. Error is estimated by comparing to fine RK4 simulation with $\delta ={10}^{-4}$.

Figure 6

Poincaré sections of $\overline{H}$ with $H(Q,P)=(Q^2+1)(P^2+1)/2$ and different $\omega$ values at the surface $x=0$ on constant energy shell $\overline{H}=10$. Initial conditions are uniformly sampled in the $q,p$ plane. Note curves do intersect because each $q,p$ pair corresponds to two $y$ values.