Local structure-preserving algorithms for general multi-symplectic Hamiltonian PDEs

Abstract

Many PDEs can be recast into the general multi-symplectic formulation possessing three local conservation laws. We devote the present paper to some systematic methods, which hold the discrete versions of the local conservation laws respectively, for the general multi-symplectic PDEs. For the original problem subjected to appropriate boundary conditions, the proposed methods are globally conservative. The proposed methods are successfully applied to many one-dimensional and multi-dimensional Hamiltonian PDEs, such as KdV equation, G–P equation, Maxwell’s equations and so on. Numerical experiments are carried out to verify the theoretical analysis.

Introduction

A wide range of $d$-dimensional PDEs can be written into a multi-symplectic Hamiltonian form [1]

$$M{\partial }_{t}z+\sum _{l=1}^d{K}_l{\partial }_{x_l}z={\nabla }_{z}S\left(z\right),z\in {\mathbb{R}}^{d_1},$$

where $M$ and $K_l$, $l=1,2,\dots ,d$, are constant $d_1\times d_1$ skew-symmetric matrices, $z=z(t,\mathbf{x})$, $\mathbf{x}=(x_1,x_2,\dots ,x_d)$, is a state variable, and $S:{\mathbb{R}}^{d_1}\to \mathbb{R}$ is a scalar-valued smooth function, such as the Sine-Gordon equation [2], the Korteweg–de Vries (KdV) equation [3], the Camassa–Holm equation [4], the regularized long-wave (RLW) equation [5], the nonlinear Schrödinger (NLS) equation [[6], [7]], Maxwell’s equations [[8], [9]] and so on. It is well-known that the solution $z$ of (1) admits a multi-symplectic conservation law (MSCL)

$${\partial }_t\omega +\sum _{l=1}^d{\partial }_{x_l}{\kappa }_l=0,$$

with the differential 2-form $\omega =dz\wedge \widehat{M}dz$ and ${\kappa }_l=dz\wedge {\widehat{K}}_{l}dz$ where matrices $\widehat{M}$ and ${\widehat{K}}_l$ satisfy

$$M=\widehat{M}-{\widehat{M}}^T,K_l={\widehat{K}}_l-{\widehat{K}}_l^T.$$

A numerical method that satisfies a discrete version of (2) is called multi-symplectic integrator; see [2] for reviews of multi-symplectic integration. Up to now, some multi-symplectic integrators have been proposed for the multi-symplectic Hamiltonian PDEs and exhibit good numerical performances, such as the box scheme obtained by applying simply implicit midpoint rule (a Runge–Kutta method) in both time and space [2], the Euler-box scheme developed by using the non-compact Euler rule in both time and space [10], the pseudospectral scheme obtained by applying the pseudospectral method in space and the implicit midpoint rule in time [11], and the diamond scheme [12]. Although multi-symplecticity by itself does not ensure good performance on traditional criteria like accuracy and stability, multi-symplectic methods do have a number of advantages. On the one hand, they are essentially variational and the standard discrete model in the relevant part of physics. On the other hand, the symplecticity in space is essential for a semidiscretization to be amenable to symplectic time integration, whose advantages in long-time integration are well known, and the symplecticity in space can also preserve periodic, quasi-periodic and heteroclinic solutions. However, multi-symplectic integrators can only preserve some conservation laws, and these do not typically include those of energy and momentum, which can be important for nonlinear stability and the nonlinear convergence analysis of the schemes.

In recent years, there has been an increased emphasis on constructing numerical methods preserving certain invariant quantities such as energy and momentum in the continuous dynamical systems. Nowadays, for computing physical/chemical problem, whether the numerical scheme captures the energy/momentum of the system has been a criterion to judge the success of the numerical simulation. Many PDEs with appropriate boundary conditions admit the global energy/momentum conservation law and there have been some techniques for constructing the global energy-/momentum-preserving methods for these conservative PDEs [[13], [14], [15]]. However, these discrete preserving properties are global and depend on the boundary conditions inevitably, so that the schemes are invalid for the problem without appropriate boundary conditions. As we know, besides possessing the global conservation laws, many PDEs also have local conservation laws which are independent of the boundary conditions and more essential than the global ones. For example, the multi-symplectic Hamiltonian PDE admits a local energy conservation law (LECL)

$${\partial }_{t}E+\sum _{l=1}^d{\partial }_{x_l}{F}_l=0,$$

where $E=S\left(z\right)+{\sum }_{l=1}^d{\left({\partial }_{x_l}z\right)}^T{\widehat{K}}_{l}z$ and $F_l=-{\left({\partial }_{t}z\right)}^T{\widehat{K}}_{l}z$ represent the energy density and energy fluxes respectively, and also has local momentum conservation laws (LMCLs) in the $x_l$ directions , $l=1,2,\dots ,d$, given by

$${\partial }_t{I}_l+{\partial }_{x_l}{G}_l+\sum _{j=1,j\ne l}^d{\partial }_{x_j}{\tilde{G}}_{l,j}=0,$$

where $I_l=-{\left({\partial }_{x_l}z\right)}^T\widehat{M}z$, $G_l=S\left(z\right)+{\left({\partial }_{t}z\right)}^T\widehat{M}z+{\sum }_{j=1,j\ne l}^d{\left({\partial }_{x_j}z\right)}^T{\widehat{K}}_{j}z$ and ${\tilde{G}}_{l,j}=-{\left({\partial }_{x_l}z\right)}^T{\widehat{K}}_{j}z$. It is clear that the local conservation laws (4), (5) are established on any time–space point and independent of boundary conditions. Assume that the PDEs are imposed on some appropriate boundary conditions, these LECL and LMCLs will result in global energy and momentum conservation laws, respectively. And thus the LECL and LMCLs produce richer informations of the PDE system than the corresponding global ones. Naturally, we expect to propose a scheme that satisfies a discrete version of LECL (4) or LMCL (5). In this paper, this kind of schemes will be called local energy-preserving (LEP) integrator or local momentum-preserving (LMP) integrator. On the construction of the LEM and LMP integrators for the PDEs, we have done some works. In [[16], [17]], by using the concatenating method, we have proposed some LEP and LMP schemes for the Sine-Gordon equation and the coupled NLS equations. These schemes have excellent properties and exhibit good numerical performance, but they are not completely systematic either in their derivations or in their applicability to a class of PDEs. Recently, combining the implicit midpoint rule with the averaged vector field (AVF) method, we have developed a LEP integrator and a LMP integrator which are applicable to the entire class (1) in one dimension [18]. The derivation of the LEP and LMP schemes are easy, however, these schemes have some less positive features. The schemes are fully implicit, which makes them expensive; The implicit equations may not have a solution: with periodic boundary conditions, solvability requires that the number of grid points be odd [19]; For a given PDE, it is difficult to establish the error estimates of the obtained LEP/LMP scheme no matter whether there exist positive invariants or not. The reason is that there is no relationship $c_1\Vert e^n\Vert \le \Vert A_x{e}^n\Vert \le c_2\Vert e^n\Vert$, where constants $c_1$, $c_2>0$, $\Vert \cdot \Vert$ represents the discrete $l_2$-norm, $A_x$ is an average operator, and $e^n$ represents the error vector whose entries denote the difference between the numerical solution of the LEP/LMP scheme and the exact solution of original problem.

In this paper, instead of discretizing one particular PDE, we wish to develop methods that are applicable to the entire class of multi-symplectic Hamiltonian PDEs (1), specializing to a particular equation or family as late as possible. To this end, we split the skew-symmetric matrix $M$ or $K_l$ in the form of (3). And then applying the leap-frog rule to space/time gives a semi-discrete system of ODEs with respect to time/space. Next, we use a discrete gradient (DG) method such as the coordinate increment DG method [20], the midpoint DG method [21], or the AVF method to integrate the obtained ODEs. In this study, an application of the AVF method gives a fully discrete LEP/LMP scheme for the multi-symplectic Hamiltonian PDE (1). The AVF method was first written down in [22] for the system of ODEs $dy∕dt=f\left(y\right)$, $y\in {\mathbb{R}}^d$, and identified as an energy-preserving and B-series method in [23]. The second-order AVF method is the map $y_n⟼y_{n+1}$ defined by

$$\frac{y_{n+1}-y_n}{\tau }={\int }_0^{1}f((1-\xi )y_n+\xi y_{n+1})\text{d}\xi ,$$

where $\tau$ is the step size. The method (6) is affine-covariant [24] and self-adjoint. When $f$ is Hamiltonian with respect to a constant symplectic structure, i.e., $f=S^{-1}\nabla H$ with $S$ a nonsingular and antisymmetric matrix, the AVF method preserves the Hamiltonian $H:{\mathbb{R}}^{d_1}\to \mathbb{R}$. Besides preserving the LECL or LMCL of the system (1), the present method also has several appealing properties. It is systematic and applicable to a huge of PDEs which has the multi-symplectic Hamiltonian form (1); Applications of the matrix splitting technique and leap-frog rule may result in a fully explicit or linearly implicit scheme for some nonlinear PDEs; Different schemes can be proposed due to different splitting matrices; The method is effective and valid no matter whether the number of the grid points is odd or even; The present method is established on the grid points unlike the method in [18] establishing on the central box, so that the discrete conservation laws are also related to the grid points. For some particular multi-symplectic Hamiltonian PDEs holding the positive discrete conservation laws, the feature allows us to establish the error estimate of the scheme easily.

The present work is organized as follows. In Section 2, we briefly review the multi-symplectic form of the one-dimensional PDE and then develop a class of energy- and momentum-preserving integrators for it. We apply the developed local structure-preserving integrators to some classical one-dimensional PDEs and conduct numerical experiments to show their numerical performance in the same section. In Section 3, the methodology of the localstructure-preserving algorithms for the one-dimensional multi-symplectic Hamiltonian PDE is extended to multi-dimensional multi-symplectic Hamiltonian PDEs. In Section 4, we study the LEP discretization by means of numerical dispersion relation for the linear PDE. Finally, we finish the paper with concluding remarks in Section 5.