Multisymplectic integrator

A partial differential equation (PDE) is said to be a multisymplectic equation if it can be written in the form $${\displaystyle Kz_{t}+Lz_{x}=\nabla S(z),}$$ where ${\displaystyle z(t,x)}$ is the unknown, ${\displaystyle K}$ and ${\displaystyle L}$ are (constant) skew-symmetric matrices and ${\displaystyle \nabla S}$ denotes the gradient of ${\displaystyle S}$.^[1] This is a natural generalization of ${\displaystyle Jz_{t}=\nabla H(z)}$, the form of a Hamiltonian ODE.^[2]

Examples of multisymplectic PDEs include the nonlinear Klein–Gordon equation ${\displaystyle u_{tt}-u_{xx}=V'(u)}$, or more generally the nonlinear wave equation ${\displaystyle u_{tt}=\partial _{x}\sigma '(u_{x})-f'(u)}$,^[3] and the KdV equation ${\displaystyle u_{t}+uu_{x}+u_{xxx}=0}$.^[4]

Define the 2-forms ${\displaystyle \omega }$ and ${\displaystyle \kappa }$ by $${\displaystyle \omega (u,v)=\langle Ku,v\rangle \quad {\text{and}}\quad \kappa (u,v)=\langle Lu,v\rangle }$$ where ${\displaystyle \langle \,\cdot \,,\,\cdot \,\rangle }$ denotes the dot product. The differential equation preserves symplecticity in the sense that^[5][6] $${\displaystyle \partial _{t}\omega +\partial _{x}\kappa =0.}$$ Taking the dot product of the PDE with ${\displaystyle u_{t}}$ yields the local conservation law for energy:^[7] $${\displaystyle \partial _{t}E(u)+\partial _{x}F(u)=0}$$ where $${\displaystyle {\begin{aligned}E(u)&=S(u)-{\tfrac {1}{2}}\kappa (u_{x},u),\\[1ex]F(u)&={\tfrac {1}{2}}\kappa (u_{t},u).\end{aligned}}}$$ The local conservation law for momentum is derived similarly:^[7] $${\displaystyle \partial _{t}I(u)+\partial _{x}G(u)=0}$$ where $${\displaystyle {\begin{aligned}I(u)&={\tfrac {1}{2}}\omega (u_{x},u),\\[1ex]G(u)&=S(u)-{\tfrac {1}{2}}\omega (u_{t},u).\end{aligned}}}$$

A multisymplectic integrator is a numerical method for solving multisymplectic PDEs whose numerical solution conserves a discrete form of symplecticity.^[8] One example is the Euler box scheme, which is derived by applying the symplectic Euler method to each independent variable.^[9]

The Euler box scheme uses a splitting of the skew-symmetric matrices ${\displaystyle K}$ and ${\displaystyle L}$ of the form: $${\displaystyle {\begin{aligned}K&=K_{+}+K_{-}&{\text{with}}&&K_{-}&=-K_{+}^{T},\\L&=L_{+}+L_{-}&{\text{with}}&&L_{-}&=-L_{+}^{T}.\end{aligned}}}$$ For instance, one can take ${\displaystyle K_{+}}$ and ${\displaystyle L_{+}}$ to be the upper triangular part of ${\displaystyle K}$ and ${\displaystyle L}$, respectively.^[10]

Now introduce a uniform grid and let ${\displaystyle u_{n,i}}$ denote the approximation to ${\displaystyle u(n\Delta {t},i\Delta {x})}$ where ${\displaystyle \Delta {t}}$ and ${\displaystyle \Delta {x}}$ are the grid spacing in the time- and space-direction. Then the Euler box scheme is $${\displaystyle K_{+}\partial _{t}^{+}u_{n,i}+K_{-}\partial _{t}^{-}u_{n,i}+L_{+}\partial _{x}^{+}u_{n,i}+L_{-}\partial _{x}^{-}u_{n,i}=\nabla {S}(u_{n,i})}$$ where the finite difference operators are defined by^[11] $${\displaystyle {\begin{aligned}\partial _{t}^{+}u_{n,i}&={\frac {u_{n+1,i}-u_{n,i}}{\Delta {t}}},&\partial _{x}^{+}u_{n,i}&={\frac {u_{n,i+1}-u_{n,i}}{\Delta {x}}},\\[1ex]\partial _{t}^{-}u_{n,i}&={\frac {u_{n,i}-u_{n-1,i}}{\Delta {t}}},&\partial _{x}^{-}u_{n,i}&={\frac {u_{n,i}-u_{n,i-1}}{\Delta {x}}}.\end{aligned}}}$$ The Euler box scheme is a first-order method,^[9] which satisfies the discrete conservation law^[12] $${\displaystyle \partial _{t}^{+}\omega _{n,i}+\partial _{x}^{+}\kappa _{n,i}=0}$$ where $${\displaystyle {\begin{aligned}\omega _{n,i}&=\mathrm {d} u_{n,i-1}\wedge K_{+}\,\mathrm {d} u_{n,i}\\[1ex]\kappa _{n,i}&=\mathrm {d} u_{n-1,i}\wedge L_{+}\,\mathrm {d} u_{n,i}.\end{aligned}}}$$