Symplectic integrator

In mathematics, a symplectic integrator (SI) is a numerical integration scheme for Hamiltonian systems. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations. They are widely used in nonlinear dynamics, molecular dynamics, discrete element methods, accelerator physics, plasma physics, quantum physics, and celestial mechanics.

Introduction

Symplectic integrators are designed for the numerical solution of Hamilton's equations, which read

{\dot {p}}=-{\frac {\partial H}{\partial q}}\quad {\mbox{and}}\quad {\dot {q}}={\frac {\partial H}{\partial p}},

where $q$ denotes the position coordinates, $p$ the momentum coordinates, and $H$ is the Hamiltonian. The set of position and momentum coordinates $(q,p)$ are called canonical coordinates. (See Hamiltonian mechanics for more background.)

The time evolution of Hamilton's equations is a symplectomorphism, meaning that it conserves the symplectic two-form $dp\wedge dq$ . A numerical scheme is a symplectic integrator if it also conserves this two-form.

Symplectic integrators possess, as a conserved quantity, a Hamiltonian which is slightly perturbed from the original one. By virtue of these advantages, the SI scheme has been widely applied to the calculations of long-term evolution of chaotic Hamiltonian systems ranging from the Kepler problem to the classical and semi-classical simulations in molecular dynamics.

Most of the usual numerical methods, like the primitive Euler scheme and the classical Runge–Kutta scheme, are not symplectic integrators.

Splitting methods for separable Hamiltonians

A widely used class of symplectic integrators is formed by the splitting methods.

Assume that the Hamiltonian is separable, meaning that it can be written in the form

H(p,q)=T(p)+V(q).\qquad \qquad (1)

This happens frequently in Hamiltonian mechanics, with T being the kinetic energy and V the potential energy.

For the notational simplicity, let us introduce the symbol $z=(q,p)$ to denote the canonical coordinates including both the position and momentum coordinates. Then, the set of the Hamilton's equations given in the introduction can be expressed in a single expression as

{\dot {z}}=\{z,H(z)\},\qquad \qquad (2)

where $\{\cdot ,\cdot \}$ is a Poisson bracket. Furthermore, by introducing an operator $D_{H}\cdot =\{\cdot ,H\}$ , which returns a Poisson bracket of the operand with the Hamiltonian, the expression of the Hamilton's equation can be further simplified to

{\dot {z}}=D_{H}z.

The formal solution of this set of equations is given as a matrix exponential:

z(\tau )=\exp(\tau D_{H})z(0).\qquad \qquad (3)

Note the positivity of $\tau D_{H}$ in the matrix exponential.

When the Hamiltonian has the form of eq. (1), the solution (3) is equivalent to

z(\tau )=\exp[\tau (D_{T}+D_{V})]z(0).\qquad \qquad (4)

The SI scheme approximates the time-evolution operator $\exp[\tau (D_{T}+D_{V})]$ in the formal solution (4) by a product of operators as

{\begin{aligned}\exp[\tau (D_{T}+D_{V})]&=\prod _{i=1}^{k}\exp(c_{i}\tau D_{T})\exp(d_{i}\tau D_{V})+O(\tau ^{k+1})\\&=\exp(c_{1}\tau D_{T})\exp(d_{1}\tau D_{V})\dots \exp(c_{k}\tau D_{T})\exp(d_{k}\tau D_{V})+O(\tau ^{k+1}),\end{aligned}}\qquad \qquad (5)

where $c_{i}$ and $d_{i}$ are real numbers, $k$ is an integer, which is called the order of the integrator, and where $\sum _{i=1}^{k}c_{i}=\sum _{i=1}^{k}d_{i}=1$ . Note that each of the operators $\exp(c_{i}\tau D_{T})$ and $\exp(d_{i}\tau D_{V})$ provides a symplectic map, so their product appearing in the right-hand side of (5) also constitutes a symplectic map.

Since $D_{T}^{2}z=\{\{z,T\},T\}=\{({\dot {q}},0),T\}=(0,0)$ for all $z$ , we can conclude that

D_{T}^{2}=0.\qquad \qquad (6)

By using a Taylor series, $\exp(aD_{T})$ can be expressed as

\exp(aD_{T})=\sum _{n=0}^{\infty }{\frac {(aD_{T})^{n}}{n!}},\qquad \qquad (7)

where $a$ is an arbitrary real number. Combining (6) and (7), and by using the same reasoning for $D_{V}$ as we have used for $D_{T}$ , we get

\left\{{\begin{aligned}\exp(aD_{T})&=1+aD_{T},\\\exp(aD_{V})&=1+aD_{V}.\end{aligned}}\right.\qquad \qquad (8)

In concrete terms, $\exp(c_{i}\tau D_{T})$ gives the mapping

{\begin{pmatrix}q\\p\end{pmatrix}}\mapsto {\begin{pmatrix}q+\tau c_{i}{\frac {\partial T}{\partial p}}(p)\\p\end{pmatrix}},

and $\exp(d_{i}\tau D_{V})$ gives

{\begin{pmatrix}q\\p\end{pmatrix}}\mapsto {\begin{pmatrix}q\\p-\tau d_{i}{\frac {\partial V}{\partial q}}(q)\\\end{pmatrix}}.

Note that both of these maps are practically computable.

Splitting methods for general nonseparable Hamiltonians

General nonseparable Hamiltonians can also be explicitly and symplectically integrated.

To do so, Tao introduced a restraint that binds two copies of phase space together to enable an explicit splitting of such systems.^[1] The idea is, instead of $H(Q,P)$ , one simulates ${\bar {H}}(q,p,x,y)=H(q,y)+H(x,p)+\omega \left(\|q-x\|_{2}^{2}/2+\|p-y\|_{2}^{2}/2\right)$ , whose solution agrees with that of $H(Q,P)$ in the sense that $q(t)=x(t)=Q(t),p(t)=y(t)=P(t)$ .

The new Hamiltonian is advantageous for explicit symplectic integration, because it can be split into the sum of three sub-Hamiltonians, $H_{A}=H(q,y)$ , $H_{B}=H(x,p)$ , and $H_{C}=\omega \left(\|q-x\|_{2}^{2}/2+\|p-y\|_{2}^{2}/2\right)$ . Exact solutions of all three sub-Hamiltonians can be explicitly obtained: both $H_{A},H_{B}$ solutions correspond to shifts of mismatched position and momentum, and $H_{C}$ corresponds to a linear transformation. To symplectically simulate the system, one simply composes these solution maps.

Examples

The simplified form of the equations are:

p_{i+1}=p_{i}+c_{i}F(q_{i})*t

q_{i+1}=q_{i}+d_{i}p_{i+1}*t

After converting into Lagrangian coordinates:

v_{i+1}=v_{i}+c_{i}{\tfrac {F(x_{i})}{m}}*t

x_{i+1}=x_{i}+d_{i}v_{i+1}*t

Where $F(x)$ is the force vector at $x$ , and $m$ is the scalar quantity of mass.

Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position $q$ and velocity $p$ .

To apply a timestep with values $c_{1,2,3},d_{1,2,3}$ to the particle, carry out the following steps:

Update the velocity of the particle by adding to it its acceleration multiplied by $c_{1}$
Update the position of the particle by adding to it its (updated) velocity multiplied by $d_{1}$
Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by $c_{2}$
Update the position of the particle by adding to it its (double-updated) velocity multiplied by $d_{2}$
Update the velocity of the particle by adding to it its acceleration (at the double-updated position) multiplied by $c_{3}$
Update the position of the particle by adding to it its (triple-updated) velocity multiplied by $d_{3}$

A first-order example

The symplectic Euler method is the first-order integrator with $k=1$ and coefficients

c_{1}=d_{1}=1.\ \

Note that the algorithm above does not work if time-reversibility is needed. The algorithm has to be implemented in two parts, one for positive time steps, one for negative time steps.

A second-order example

The Verlet method is the second-order integrator with $k=2$ and coefficients

c_{1}=0,\qquad c_{2}=1,\qquad d_{1}=d_{2}={\tfrac {1}{2}}.

Since $c_{1}=0$ , the algorithm above is symmetric in time. There is 3 steps to the algorithm, and step 1 and 3 are exactly the same, so the positive time version can be used for negative time.

A third-order example

A third-order symplectic integrator (with $k=3$ ) was discovered by Ronald Ruth in 1983.^[2] One of the many solutions is given by

{\begin{aligned}c_{1}&=1,&c_{2}&=-{\tfrac {2}{3}},&c_{3}&={\tfrac {2}{3}},\\d_{1}&=-{\tfrac {1}{24}},&d_{2}&={\tfrac {3}{4}},&d_{3}&={\tfrac {7}{24}}.\end{aligned}}

A fourth-order example

A fourth-order integrator (with $k=4$ ) was also discovered by Ruth in 1983 and distributed privately to the particle-accelerator community at that time. This was described in a lively review article by Forest.^[3] This fourth-order integrator was published in 1990 by Forest and Ruth and also independently discovered by two other groups around that same time.^[4]^[5]^[6]

{\begin{aligned}c_{1}&=c_{4}={\frac {1}{2(2-2^{1/3})}},&c_{2}&=c_{3}={\frac {1-2^{1/3}}{2(2-2^{1/3})}},\\d_{1}&=d_{3}={\frac {1}{2-2^{1/3}}},&d_{2}&=-{\frac {2^{1/3}}{2-2^{1/3}}},\quad d_{4}=0.\end{aligned}}

To determine these coefficients, the Baker–Campbell–Hausdorff formula can be used. Yoshida, in particular, gives an elegant derivation of coefficients for higher-order integrators. Later on, Blanes and Moan^[7] further developed partitioned Runge–Kutta methods for the integration of systems with separable Hamiltonians with very small error constants.

References

↑ Tao, Molei (2016). "Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm and long time performance". Phys. Rev. E. 94: 043303. arXiv:1609.02212. Bibcode:2016PhRvE..94d3303T. doi:10.1103/PhysRevE.94.043303.
↑ Ruth, Ronald D. (August 1983). "A Canonical Integration Technique". Nuclear Science, IEEE Trans. on. NS-30 (4): 2669–2671. Bibcode:1983ITNS...30.2669R. doi:10.1109/TNS.1983.4332919.
↑ Forest, Etienne (2006). "Geometric Integration for Particle Accelerators". J. Phys. A: Math. Gen. 39 (19): 5321–5377. Bibcode:2006JPhA...39.5321F. doi:10.1088/0305-4470/39/19/S03.
↑ Forest, E.; Ruth, Ronald D. (1990). "Fourth-order symplectic integration". Physica D. 43: 105. Bibcode:1990PhyD...43..105F. doi:10.1016/0167-2789(90)90019-L.
↑ Yoshida, H. (1990). "Construction of higher order symplectic integrators". Phys. Lett. A. 150 (5–7): 262. Bibcode:1990PhLA..150..262Y. doi:10.1016/0375-9601(90)90092-3.
↑ Candy, J.; Rozmus, W (1991). "A Symplectic Integration Algorithm for Separable Hamiltonian Functions". J. Comput. Phys. 92: 230. Bibcode:1991JCoPh..92..230C. doi:10.1016/0021-9991(91)90299-Z.
↑ Blanes, S.; Moan, P. C. (May 2002). "Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods". Journal of Computational and Applied Mathematics. 142 (2): 313–330. Bibcode:2002JCoAM.142..313B. doi:10.1016/S0377-0427(01)00492-7.

Leimkuhler, Ben; Reich, Sebastian (2005). Simulating Hamiltonian Dynamics. Cambridge University Press. ISBN 0-521-77290-7.
Hairer, Ernst; Lubich, Christian; Wanner, Gerhard (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (2 ed.). Springer. ISBN 978-3-540-30663-4.

This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[1] Tao, Molei (2016). "Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm and long time performance". Phys. Rev. E. 94: 043303. arXiv:1609.02212. Bibcode:2016PhRvE..94d3303T. doi:10.1103/PhysRevE.94.043303.

[2] Ruth, Ronald D. (August 1983). "A Canonical Integration Technique". Nuclear Science, IEEE Trans. on. NS-30 (4): 2669–2671. Bibcode:1983ITNS...30.2669R. doi:10.1109/TNS.1983.4332919.

[3] Forest, Etienne (2006). "Geometric Integration for Particle Accelerators". J. Phys. A: Math. Gen. 39 (19): 5321–5377. Bibcode:2006JPhA...39.5321F. doi:10.1088/0305-4470/39/19/S03.

[4] Forest, E.; Ruth, Ronald D. (1990). "Fourth-order symplectic integration". Physica D. 43: 105. Bibcode:1990PhyD...43..105F. doi:10.1016/0167-2789(90)90019-L.

[5] Yoshida, H. (1990). "Construction of higher order symplectic integrators". Phys. Lett. A. 150 (5–7): 262. Bibcode:1990PhLA..150..262Y. doi:10.1016/0375-9601(90)90092-3.

[6] Candy, J.; Rozmus, W (1991). "A Symplectic Integration Algorithm for Separable Hamiltonian Functions". J. Comput. Phys. 92: 230. Bibcode:1991JCoPh..92..230C. doi:10.1016/0021-9991(91)90299-Z.

[7] Blanes, S.; Moan, P. C. (May 2002). "Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods". Journal of Computational and Applied Mathematics. 142 (2): 313–330. Bibcode:2002JCoAM.142..313B. doi:10.1016/S0377-0427(01)00492-7.