Symplectic integrator

1 Hamiltonian mechanics

1.1 Generalized Momentum

The generalized momentum $p_{i}$ corresponding to the generalized coordinate $q^{i}$ is given by:

p_{i}=\partialderivative{L}{\dot{q}^{i}}

where $L$ is the Lagrangian.

1.2 Hamiltonian

Hamiltonian $H$ is given by:

H(q,p,t)=\sum_{i}p_{i}\,\dot{q}^{i}(q,p)-L(q,\dot{q}(q,p),t)

where $L(q,\dot{q}(q,p),t)$ is the Lagrangian.

1.3 Hamilton’s Equations

	$\displaystyle\dot{q}^{i}$	$\displaystyle=\partialderivative{H}{p_{i}}$
	$\displaystyle\dot{p}_{i}$	$\displaystyle=-\partialderivative{H}{q^{i}}$

Or in vector notation,

\boldsymbol{z}=\begin{pmatrix}\boldsymbol{q}\\ \boldsymbol{p}\end{pmatrix}\,,\ J=\begin{pmatrix}0&I\\ -I&0\end{pmatrix}

\dot{\boldsymbol{z}}=J\gradient_{\boldsymbol{z}}{H(\boldsymbol{z})}

2 Poisson Bracket

\{A(q,p,t),B(q,p,t)\}=\sum_{i}\left(\partialderivative{A}{q^{i}}% \partialderivative{B}{p_{i}}-\partialderivative{B}{q^{i}}\partialderivative{A}% {p_{i}}\right)

\derivative{A(q,p,t)}{t}=\{A,H\}+\partialderivative{A}{t}

3 Symplectic Condition

Under a canonical transformation ( $\boldsymbol{z}\to\boldsymbol{w}$ ) :

	$\displaystyle\dot{\boldsymbol{z}}$	$\displaystyle=J\gradient_{\boldsymbol{z}}{H(\boldsymbol{z})}$
	$\displaystyle\dot{\boldsymbol{w}}$	$\displaystyle=J\gradient_{\boldsymbol{w}}{H(\boldsymbol{w})}$

The Jacobian matrix $M$ is defined as

M_{ij}=\partialderivative{w_{i}}{z_{j}}

Using the chain rule:

	$\displaystyle(\gradient_{\boldsymbol{z}}{H})_{i}$	$\displaystyle=\partialderivative{H}{z_{i}}$
	$\displaystyle\partialderivative{H}{z_{i}}$	$\displaystyle=\sum_{j}\partialderivative{H}{w_{j}}\partialderivative{w_{j}}{z_% {i}}$
		$\displaystyle=\sum_{j}(M^{T})_{ij}(\gradient_{w}H)_{j}$
	$\displaystyle\gradient_{\boldsymbol{z}}{H(\boldsymbol{z})}$	$\displaystyle=M^{T}\gradient_{w}{H(\boldsymbol{w})}$

	$\displaystyle\derivative{w_{i}}{t}$	$\displaystyle=\sum_{j}\partialderivative{w_{i}}{z_{j}}\derivative{z_{j}}{t}$
	$\displaystyle\dot{w_{i}}$	$\displaystyle=\sum_{j}M_{ij}\dot{z_{j}}$
	$\displaystyle\dot{\boldsymbol{w}}$	$\displaystyle=M\dot{\boldsymbol{z}}$
		$\displaystyle=MJ\gradient_{\boldsymbol{z}}{H(\boldsymbol{z})}$
		$\displaystyle=MJM^{T}\gradient_{\boldsymbol{w}}{H(\boldsymbol{w})}$

Canonical transformation $\iff$ the symplectic condition:

MJM^{T}=J

4 Splitting methods for separable Hamiltonians

\hat{D}_{f}=\{\cdot,f\}

\dot{\boldsymbol{z}}=\{\boldsymbol{z},H\}=\hat{D}_{H}\boldsymbol{z}

(1)

By Taylor expansion:

	$\displaystyle\boldsymbol{z}(t+\Delta t)$	$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\derivative[n]{\boldsymbol{z}(t)}% {t}\left(\Delta t\right)^{n}$
		$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}(\Delta t)^{n}\hat{D}_{H}^{n}\,% \boldsymbol{z}(t)$
		$\displaystyle=\exp(\Delta t\hat{D}_{H})\boldsymbol{z}(t)$

Thus, the exact solution of (1) is given by:

\boldsymbol{z}(t)=\exp(t\hat{D}_{H})\boldsymbol{z}(0)

If the Hamiltonian is separable,

H(\boldsymbol{q},\boldsymbol{p})=K(\boldsymbol{p})+U(\boldsymbol{q})

the operator splits as:

	$\displaystyle\hat{D}_{H}$	$\displaystyle=\{\cdot,H\}=\{\cdot,K+U\}$
		$\displaystyle=\{\cdot,K\}+\{\cdot,U\}$
		$\displaystyle=\hat{D}_{K}+\hat{D}_{U}$

\exp(\Delta t\hat{D}_{H})=\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))=\exp(\Delta t% \hat{D}_{K}+\Delta t\hat{D}_{U})

Because $\hat{D}_{K}$ and $\hat{D}_{U}$ are non-commutative operators ( $[\hat{D}_{K},\hat{D}_{U}]\neq 0$ ),

\exp(\Delta t\hat{D}_{K}+\Delta t\hat{D}_{U})\neq\exp(\Delta t\hat{D}_{K})\exp% (\Delta t\hat{D}_{U})

(refer to the Baker–Campbell–Hausdorff formula)

Instead, the time-evolution operator $\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))$ can be approximated by a product of operators:

\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))=\prod_{i=k,\cdots,1}\exp(d_{i}\Delta t% \hat{D}_{K})\exp(c_{i}\Delta t\hat{D}_{U})+\mathcal{O}\big{(}(\Delta t)^{k+1}% \big{)}

(2)

where $k\in\mathbb{N}$ represents the order of the integrator, and $c_{i},d_{i}\in\mathbb{R}\ (i=1,\cdots,k)$ are coefficients chosen such that the residual error is of the order of $(\Delta t)^{k+1}$

In a direct approach, conditions for the coefficients are found by comparing the Taylor expansions of both sides of (2)

	$\displaystyle\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))$	$\displaystyle=\sum_{n=0}^{k}\frac{1}{n!}(\Delta t)^{n}(\hat{D}_{K}+\hat{D}_{U}% )^{n}+\mathcal{O}\big{(}(\Delta t)^{k+1}\big{)}$
	$\displaystyle\exp(d_{i}\Delta t\hat{D}_{K})$	$\displaystyle=\sum_{n=0}^{k}\frac{1}{n!}d_{i}^{n}(\Delta t)^{n}\hat{D}_{K}^{n}% +\mathcal{O}\big{(}(\Delta t)^{k+1}\big{)}$
	$\displaystyle\exp(c_{i}\Delta t\hat{D}_{U})$	$\displaystyle=\sum_{n=0}^{k}\frac{1}{n!}c_{i}^{n}(\Delta t)^{n}\hat{D}_{U}^{n}% +\mathcal{O}\big{(}(\Delta t)^{k+1}\big{)}$

4.1 First-Order Integrator ( $k=1$ )

	$\displaystyle\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))$	$\displaystyle=1+\Delta t(\hat{D}_{K}+\hat{D}_{U})+\mathcal{O}\big{(}(\Delta t)% ^{2}\big{)}$
	$\displaystyle\exp(d_{i}\Delta t\hat{D}_{K})$	$\displaystyle=1+d_{i}\Delta t\hat{D}_{K}+\mathcal{O}\big{(}(\Delta t)^{2}\big{)}$
	$\displaystyle\exp(c_{i}\Delta t\hat{D}_{U})$	$\displaystyle=1+c_{i}\Delta t\hat{D}_{U}+\mathcal{O}\big{(}(\Delta t)^{2}\big{)}$
	$\displaystyle\exp(d_{1}\Delta t\hat{D}_{K})\exp(c_{1}\Delta t\hat{D}_{U})$	$\displaystyle=1+\Delta t(d_{1}\hat{D}_{K}+c_{1}\hat{D}_{U})+\mathcal{O}\big{(}% (\Delta t)^{2}\big{)}$

The solution is Lie product formula:

c_{1}=d_{1}=1

\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))=\exp(\Delta t\hat{D}_{K})\exp(\Delta t% \hat{D}_{U})+\mathcal{O}\big{(}(\Delta t)^{2}\big{)}

This corresponds to the Symplectic Euler method.

4.2 Second-Order Integrator ( $k=2$ )

	$\displaystyle\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))$	$\displaystyle=1+\Delta t(\hat{D}_{K}+\hat{D}_{U})+\frac{1}{2}(\Delta t)^{2}(% \hat{D}_{K}^{2}+\hat{D}_{U}\hat{D}_{K}+\hat{D}_{K}\hat{D}_{U}+\hat{D}_{U}^{2})% +\mathcal{O}\big{(}(\Delta t)^{3}\big{)}$
	$\displaystyle\exp(d_{i}\Delta t\hat{D}_{K})$	$\displaystyle=1+d_{i}\Delta t\hat{D}_{K}+\frac{1}{2}d_{i}^{2}(\Delta t)^{2}% \hat{D}_{K}^{2}+\mathcal{O}\big{(}(\Delta t)^{3}\big{)}$
	$\displaystyle\exp(c_{i}\Delta t\hat{D}_{U})$	$\displaystyle=1+c_{i}\Delta t\hat{D}_{U}+\frac{1}{2}c_{i}^{2}(\Delta t)^{2}% \hat{D}_{U}^{2}+\mathcal{O}\big{(}(\Delta t)^{3}\big{)}$
	$\displaystyle\prod_{i=2,1}\exp(d_{i}\Delta t\hat{D}_{K})\exp(c_{i}\Delta t\hat% {D}_{U})$	$\displaystyle=1+\Delta t\left((d_{1}+d_{2})\hat{D}_{K}+(c_{1}+c_{2})\hat{D}_{U% }\right)$
		$\displaystyle\quad+(\Delta t)^{2}\left(\frac{1}{2}(d_{1}+d_{2})^{2}\hat{D}_{K}% ^{2}+\frac{1}{2}(c_{1}+c_{2})^{2}\hat{D}_{U}^{2}+(d_{1}c_{1}+d_{2}c_{2})\hat{D% }_{K}\hat{D}_{U}+c_{2}d_{1}\hat{D}_{U}\hat{D}_{K}\right)$

d_{1}+d_{2}=1\,,\ c_{1}+c_{2}=1\,,\ d_{1}c_{1}+d_{2}c_{2}=\frac{1}{2}\,,\ d_{1% }c_{2}=\frac{1}{2}

	$\displaystyle d_{2}=1-d_{1}\,,\ c_{2}=1-c_{1}$
	$\displaystyle d_{1}c_{1}+(1-d_{1})(1-c_{1})=\frac{1}{2}\,,\ d_{1}(1-c_{1})=% \frac{1}{2}$
	$\displaystyle 2d_{1}c_{1}-d_{1}-c_{1}=-\frac{1}{2}\,,\ d_{1}-\frac{1}{2}=d_{1}% c_{1}$
	$\displaystyle d_{1}-c_{1}=\frac{1}{2}$

Two prominent solutions emerge (Strang Splitting):

c_{1}=\frac{1}{2}\,,\ d_{1}=1\,,\ c_{2}=\frac{1}{2}\,,\ d_{2}=0

\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))=\exp(\frac{1}{2}\Delta t\hat{D}_{U})% \exp(\Delta t\hat{D}_{K})\exp(\frac{1}{2}\Delta t\hat{D}_{U})+\mathcal{O}\big{% (}(\Delta t)^{3}\big{)}

which correspons to ’kick-drift-kick’ form of leapfrog integration, and

c_{1}=0\,,\ d_{1}=\frac{1}{2}\,,\ c_{2}=1\,,\ d_{2}=\frac{1}{2}

\exp(\Delta t(\hat{D}_{K}+\hat{D}_{U}))=\exp(\frac{1}{2}\Delta t\hat{D}_{K})% \exp(\Delta t\hat{D}_{U})\exp(\frac{1}{2}\Delta t\hat{D}_{K})+\mathcal{O}\big{% (}(\Delta t)^{3}\big{)}

which correspons to ’drift-kick-drift’ form of leapfrog integration.

In general, two conditions emerge from the coefficients of $\hat{D}_{K}$ and $\hat{D}_{U}$ ,

\sum_{i=1}^{k}c_{i}=1\,,\ \sum_{i=1}^{k}d_{i}=1

and higher order terms provide additional conditions.

However, as noted by Yoshida: “with this direct method, it is almost hopeless to obtain a much higher integrator” [1]. Instead, coefficient sets for higher-order integrators (e.g., $k=4$ ) are systematically derived using the Baker–Campbell–Hausdorff formula (see Yoshida [1]).

4.3 Action on Phase-Space Vectors

Because the Hamiltonian is separable,

	$\displaystyle\dot{q}^{i}=\partialderivative{H}{p_{i}}=\partialderivative{K}{p_% {i}}$
	$\displaystyle\dot{p}_{i}=-\partialderivative{H}{q^{i}}=-\partialderivative{U}{% q^{i}}$

\displaystyle\hat{D}_{K}\boldsymbol{z}=\{\boldsymbol{z},K\}=\begin{pmatrix}\{% \boldsymbol{q},K\}\\ \{\boldsymbol{p},K\}\end{pmatrix}

	$\displaystyle\{\boldsymbol{q},K\}$	$\displaystyle=\begin{pmatrix}\{q_{1},K\}\\ \vdots\\ \{q_{n},K\}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\sum_{i}\left(\partialderivative{q_{1}}{q^{i}}% \partialderivative{K}{p_{i}}-\partialderivative{K}{q^{i}}\partialderivative{q_% {1}}{p_{i}}\right)\\ \vdots\\ \sum_{i}\left(\partialderivative{q_{n}}{q^{i}}\partialderivative{K}{p_{i}}-% \partialderivative{K}{q^{i}}\partialderivative{q_{n}}{p_{i}}\right)\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\partialderivative{K}{p_{1}}\\ \vdots\\ \partialderivative{K}{p_{n}}\end{pmatrix}=\begin{pmatrix}\dot{q}^{1}\\ \vdots\\ \dot{q}^{n}\end{pmatrix}=\dot{\boldsymbol{q}}$

	$\displaystyle\{\boldsymbol{p},K\}$	$\displaystyle=\begin{pmatrix}\{p_{1},K\}\\ \vdots\\ \{p_{n},K\}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\sum_{i}\left(\partialderivative{p_{1}}{q^{i}}% \partialderivative{K}{p_{i}}-\partialderivative{K}{q^{i}}\partialderivative{p_% {1}}{p_{i}}\right)\\ \vdots\\ \sum_{i}\left(\partialderivative{p_{n}}{q^{i}}\partialderivative{K}{p_{i}}-% \partialderivative{K}{q^{i}}\partialderivative{p_{n}}{p_{i}}\right)\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}0\\ \vdots\\ 0\end{pmatrix}=\boldsymbol{0}$

\hat{D}_{K}\boldsymbol{z}=\{\boldsymbol{z},K\}=\begin{pmatrix}\dot{\boldsymbol% {q}}\\ \boldsymbol{0}\end{pmatrix}

(3)

\hat{D}_{K}^{2}\boldsymbol{z}=\{\{\boldsymbol{z},K\},K\}=\begin{pmatrix}\{\dot% {\boldsymbol{q}},K\}\\ \{\boldsymbol{0},K\}\end{pmatrix}=\boldsymbol{0}

(4)

\displaystyle\hat{D}_{U}\boldsymbol{z}=\{\boldsymbol{z},U\}=\begin{pmatrix}\{% \boldsymbol{q},U\}\\ \{\boldsymbol{p},U\}\end{pmatrix}

	$\displaystyle\{\boldsymbol{q},U\}$	$\displaystyle=\begin{pmatrix}\{q_{1},U\}\\ \vdots\\ \{q_{n},U\}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\sum_{i}\left(\partialderivative{q_{1}}{q^{i}}% \partialderivative{U}{p_{i}}-\partialderivative{U}{q^{i}}\partialderivative{q_% {1}}{p_{i}}\right)\\ \vdots\\ \sum_{i}\left(\partialderivative{q_{n}}{q^{i}}\partialderivative{U}{p_{i}}-% \partialderivative{U}{q^{i}}\partialderivative{q_{n}}{p_{i}}\right)\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}0\\ \vdots\\ 0\end{pmatrix}=\boldsymbol{0}$

	$\displaystyle\{\boldsymbol{p},U\}$	$\displaystyle=\begin{pmatrix}\{p_{1},U\}\\ \vdots\\ \{p_{n},U\}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\sum_{i}\left(\partialderivative{p_{1}}{q^{i}}% \partialderivative{U}{p_{i}}-\partialderivative{U}{q^{i}}\partialderivative{p_% {1}}{p_{i}}\right)\\ \vdots\\ \sum_{i}\left(\partialderivative{p_{n}}{q^{i}}\partialderivative{U}{p_{i}}-% \partialderivative{U}{q^{i}}\partialderivative{p_{n}}{p_{i}}\right)\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}-\partialderivative{U}{q^{1}}\\ \vdots\\ -\partialderivative{U}{q^{n}}\end{pmatrix}=\begin{pmatrix}\dot{p}_{1}\\ \vdots\\ \dot{p}_{n}\end{pmatrix}=\dot{\boldsymbol{p}}$

\hat{D}_{U}\boldsymbol{z}=\{\boldsymbol{z},U\}=\begin{pmatrix}\boldsymbol{0}\\ \dot{\boldsymbol{p}}\end{pmatrix}

(5)

\hat{D}_{U}^{2}\boldsymbol{z}=\{\{\boldsymbol{z},U\},U\}=\begin{pmatrix}\{% \boldsymbol{0},U\}\\ \{\dot{\boldsymbol{p}},U\}\end{pmatrix}=\boldsymbol{0}

(6)

Using a Taylor expansion, for $\alpha\in\mathbb{R}$ ,

	$\displaystyle\exp(\alpha\hat{D}_{U})\boldsymbol{z}$	$\displaystyle=\sum_{n=0}^{\infty}\frac{\alpha^{n}\hat{D}_{U}^{n}}{n!}% \boldsymbol{z}=(1+\alpha\hat{D}_{U})\boldsymbol{z}$
	$\displaystyle\exp(\alpha\hat{D}_{K})\boldsymbol{z}$	$\displaystyle=\sum_{n=0}^{\infty}\frac{\alpha^{n}\hat{D}_{K}^{n}}{n!}% \boldsymbol{z}=(1+\alpha\hat{D}_{K})\boldsymbol{z}$

	$\displaystyle\exp(c_{i}\Delta t\hat{D}_{U})\begin{pmatrix}\boldsymbol{q}\\ \boldsymbol{p}\end{pmatrix}$	$\displaystyle=\begin{pmatrix}\boldsymbol{q}\\ \boldsymbol{p}+c_{i}\Delta t\dot{\boldsymbol{p}}\end{pmatrix}$
	$\displaystyle\exp(d_{i}\Delta t\hat{D}_{K})\begin{pmatrix}\boldsymbol{q}\\ \boldsymbol{p}\end{pmatrix}$	$\displaystyle=\begin{pmatrix}\boldsymbol{q}+d_{i}\Delta t\dot{\boldsymbol{q}}% \\ \boldsymbol{p}\end{pmatrix}$

References

[1] H. Yoshida (1990) Construction of higher order symplectic integrators. Physics Letters A 150 (5), pp. 262–268. External Links: ISSN 0375-9601, Document, Link Cited by: §4.2.