Linearly implicit local and global energy-preserving methods for PDEs with a cubic Hamiltonian

TL;DR

This paper introduces linearly implicit numerical methods that conserve energy for multi-symplectic PDEs with cubic invariants, demonstrating improved stability and efficiency over fully implicit schemes through tests on Korteweg-de Vries and Zakharov-Kuznetsov equations.

Contribution

The paper develops novel linearly implicit energy-preserving methods for multi-symplectic PDEs with cubic invariants, enhancing computational efficiency and stability.

Findings

Methods successfully conserve discrete energy laws.

Numerical tests show good stability and efficiency.

Methods outperform fully implicit schemes in speed.

Abstract

We present linearly implicit methods that preserve discrete approximations to local and global energy conservation laws for multi-symplectic PDEs with cubic invariants. The methods are tested on the one-dimensional Korteweg-de Vries equation and the two-dimensional Zakharov-Kuznetsov equation; the numerical simulations confirm the conservative properties of the methods, and demonstrate their good stability properties and superior running speed when compared to fully implicit schemes.

Tables

Table 1. Table 1 : Computational time, in seconds, for finding the solution of the first test problem at time t = 5 𝑡 5 t=5 by a temporal step size Δ ​ t = 0.005 Δ 𝑡 0.005 \Delta t=0.005 and various number of discretization points in space, M 𝑀 M .

$M$	$200$	$400$	$600$	$800$	$1000$	$1500$	$2000$
LEP	1.87	3.16	4.43	11.18	13.81	21.53	28.54
LILEP	4.24e-1	7.40e-1	1.07	1.39	1.73	2.67	3.58
GEP	12.29	78.11	242.48	1016.57	1888.69	5793.18	13154.20
LIGEP	2.16	11.15	33.50	73.94	136.93	398.53	894.52

Table 2. Table 2 : Phase and shape errors and the computational time (CT) for different schemes applied to test problem 2 of the KdV equation, for varying number of discretization points M 𝑀 M , with time step Δ ​ t = 0.01 Δ 𝑡 0.01 \Delta t=0.01 and end time t = 100 𝑡 100 t=100 .

$M$	200			400			600
	${ϵ}_{\text{shape}}$	${ϵ}_{\text{phase}}$	CT	${ϵ}_{\text{shape}}$	${ϵ}_{\text{phase}}$	CT	${ϵ}_{\text{shape}}$	${ϵ}_{\text{phase}}$	CT
LEP	4.67e-3	1.12	21.86	1.22e-3	3.81e-1	35.89	5.86e-4	2.43e-1	51.92
LILEP	4.10e-3	1.23	5.14	5.26e-4	4.88e-1	8.26	1.45e-4	3.50e-1	10.89
GEP	1.62e-2	8.61e-1	19.53	3.66e-3	1.16e-1	34.09	1.71e-3	2.32e-2	49.45
LIGEP	1.71e-2	7.50e-1	6.84	4.39e-3	5.19e-5	8.10	2.47e-3	1.31e-1	12.52

Table 3. Table 3 : Running time, in seconds, for computing 100 100 100 steps in time by the various schemes and various number of discretization points M = M x = M y 𝑀 subscript 𝑀 𝑥 subscript 𝑀 𝑦 M=M_{x}=M_{y} in each spatial direction, solving our test problem for the Zakharov–Kuznetsov equation.

$M$	$45$	$75$	$105$	$135$	$165$	$195$	$225$	$255$
LEP	5.10	32.20	48.43	101.59	125.23	258.64	353.98	510.00
LILEP	2.04	8.87	14.57	31.02	37.25	78.98	108.02	157.91
GEP	3.62	19.54	41.87	73.59	122.31	186.74	258.19	352.36
LIGEP	1.38	6.00	13.45	23.79	39.31	60.27	83.32	113.13

Equations

$$\dot{y}=f(y)=\hat{Q}(y)+\hat{B}y+\hat{c},y\in {\mathbb{R}}^M,$$

$$\frac{y^{n+1}-y^n}{\Delta t}=\bar{Q}(y^n,y^{n+1})+\hat{B}\frac{y^n+y^{n+1}}{2}+\hat{c},$$

$$\displaystyle\bar{Q}(y^{n},y^{n+1})=\frac{1}{2}\big{(}\hat{Q}(y^{n}+y^{n+1})-\hat{Q}(y^{n})-\hat{Q}(y^{n+1})\big{)}$$

$$\dot{y}=A\nabla H(y),$$

$$\frac{y^{n+1}-y^n}{\Delta t}=3A\bar{H}(y^n,y^{n+1},\cdot ),$$

$$\displaystyle\frac{y^{n+1}-y^{n}}{\Delta t}=3A\frac{\partial\bar{H}}{\partial x}\bigg{\lvert}_{(y^{n},y^{n+1})},$$

$$H(y)=y^{T}Q(y)y+y^{T}By+c^{T}y+d,$$

$$\dot{\tilde{y}}=\tilde{A}\nabla \tilde{H}(\tilde{y})$$

$$\bar{H}(x,y,z)=x^{T}Q(y)z+\frac{1}{3}(x^{T}By+y^{T}Bz+z^{T}Bx)+\frac{1}{3}{c}^T(x+y+z)+d.$$

$$K{z}_t+L{z}_x=\nabla S(z),z\in {\mathbb{R}}^l,(x,t)\in \mathbb{R}\times \mathbb{R},$$

$${\partial }_t\omega +{\partial }_x\kappa =0,\omega =dz\wedge K_{+}dz,\kappa =dz\wedge L_{+}dz,$$

$$E_t+F_x=0,E=S(z)+z_x^T{L}_{+}z,F=-z_t^T{L}_{+}z,$$

$$I_t+G_x=0,G=S(z)+z_t^T{K}_{+}z,I=-z_x^T{K}_{+}z,$$

$$K=K_{+}-K_{+}^T,L=L_{+}-L_{+}^T.$$

$$K{z}_t+\sum _{\alpha =1}^d{L}^{\alpha }{z}_{x_{\alpha }}=\nabla S(z),z\in {\mathbb{R}}^l,(x,t)\in {\mathbb{R}}^d\times \mathbb{R},$$

$$E_t+\sum _{\alpha =1}^d{F}_{x_{\alpha }}^{\alpha }=0,$$

$$\frac{d}{dt}\mathcal{E}(z)=0,$$

$$u_t+\eta u{u}_x+{\gamma }^2{u}_{xxx}=0,$$

$$\begin{array}{rl}\frac{1}{2}{u}_t+w_x& =0,\\ -\frac{1}{2}{\varphi }_t-\gamma v_x& =-w+\frac{\eta }{2}{u}^2,\\ \gamma u_x& =v,\\ -{\varphi }_x& =-u,\end{array}$$

$$K=\left[\begin{array}{cccc}0& \frac{1}{2}& 0& 0\\ -\frac{1}{2}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],L=\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& -\gamma & 0\\ 0& \gamma & 0& 0\\ -1& 0& 0& 0\end{array}\right].$$

$$u_t+u{u}_x+u_{xxx}+u_{xyy}=0.$$

$$K{z}_t+L^1{z}_x+L^2{z}_y=\nabla S(z),z\in {\mathbb{R}}^6,(x,y,t)\in {\mathbb{R}}^2\times \mathbb{R}.$$

$$\begin{array}{rl}{\varphi }_x& =u,\\ \frac{1}{2}{\varphi }_t+v_x+w_y& =p-\frac{1}{2}{u}^2,\\ w_x-v_y& =0,\\ -\frac{1}{2}{u}_t-p_x& =0,\\ -u_x+q_y& =-v,\\ -q_x-u_y& =-w,\end{array}$$

$$k_{2,4}=\frac{1}{2},k_{4,2}=-\frac{1}{2},$$

$$l_{1,4}^1=l_{2,5}^1=l_{3,6}^1=1,l_{4,1}^1=l_{5,2}^1=l_{6,2}^1=-1,$$

$$l_{2,6}^2=l_{5,3}^2=1,l_{6,2}^2=l_{3,5}^2=-1.$$

$${\delta }_t{v}_j^n$$

$${\mu }_t{v}_j^n$$

$${\delta }_t(uv{)}_j^n=(\epsilon u_j^{n+1}+(1-\epsilon )u_j^n){\delta }_t{v}_j^n+{\delta }_t{u}_j^n((1-\epsilon )v_j^{n+1}+\epsilon v_j^n),0\le \epsilon \le 1.$$

$${\delta }_t(uv{)}_j^n$$

$$K{\partial }_t{\mu }_x{z}_j+L{\delta }_x{z}_j=\nabla S({\mu }_x{z}_j),j=0,\dots ,M-1.$$

Full text

Keywords: Structure-preserving methods, multi-symplectic PDEs, Kahan’s method.

Classification: 37K05, 65M06, 65P10

1 Introduction

In recent years, much attention has been given to the design and analysis of numerical methods for differential equations that can capture geometric properties of the exact flow. The increased interest in this subject can mainly be attributed to the superior qualitative behaviour over long time integration of such structure-preserving methods, see [1, 2, 3]. A popular class of structure-preserving methods are energy-preserving methods. Energy preservation has a far-reaching importance throughout the physical sciences [4, 5]. In particular, it has been found to be crucial in the proof of stability for several numerical methods, see e.g. [6].

Energy-preserving methods are well studied for finite-dimensional Hamiltonian systems [7, 8, 9, 10]. It is also highly conceivable that the ideas behind the finite-dimensional setting can be extended to the infinite-dimensional Hamiltonian systems or Hamiltonian partial differential equations (PDEs) [11]. There are two popular ways to construct energy-preserving methods for Hamiltonian PDEs. One approach is to semi-discretize the PDE in space so that one obtains a system of Hamiltonian ordinary differential equations (ODEs), and then apply an energy-preserving method to this semi-discrete system, see for example [10]. In this way, it is straightforward to generalise the energy-preserving methods for finite-dimensional Hamiltonian systems to Hamiltonian PDEs. However, such methods conserve only a global energy that relies on a proper boundary condition, such as a periodic boundary condition. If this is not present, the energy-preserving property will be destroyed. The other approach is based on a reformulation of the Hamiltonian PDE into a multi-symplectic form, which provides the PDE with three local conservation laws: the multi-symplectic conservation law, the energy conservation law and the momentum conservation law [5, 12, 13]. Then one may consider methods that preserve the local conservation laws, see for example multi-symplectic integrators [14] and integrators which preserve the energy conservation law or the momentum conservation law [15, 16]. Although multi-symplectic integrators, like symplectic integrators, are proven to give good numerical approximations of the solution, such methods can only preserve quadratic invariants [17, 15]. The focus of this paper is on integrators that preserve the local energy conservation law [18]. These locally defined properties are not dependent on the choice of boundary conditions, giving the methods that preserve local energy an advantage over methods that preserve a global energy, especially since local conservation laws will always lead to global conservation laws whenever periodic boundary conditions are considered. The concept of a multi-symplectic structure for PDEs was introduced by Bridges in [5, 12], see also [19] for a framework based on a Lagrangian formulation of the Cartan form. Local energy-preserving methods were first studied in [20], and have garnered much interest recently, see for example [18, 21, 22].

Most of the local energy-preserving methods proposed so far are fully implicit methods, for which a non-linear system must be solved at each time step. This is normally done by using an iterative solver where a linear system is solved at each iteration, which can lead to computationally expensive procedures, especially since the number of iterations needed in general increases with the size of the system. A fully explicit method on the other hand, may over-simplify the problem and often has inferior stability properties, so that a strong restriction on the grid ratio is needed. A good alternative may therefore be to develop linearly implicit schemes, where the solution at the next time step is found by solving only one linear system.

One example of linearly implicit methods for Hamiltonian ODEs is Kahan’s method, which was designed for solving quadratic ODEs [23] and whose geometric properties have been studied in a series of papers by Celledoni et al. [24, 25, 26]. For Hamiltonian PDEs, Matsuo and Furihata proposed the idea of using multiple points to discretize the variational derivative and thus design linearly implicit energy-preserving schemes [27]. Dahlby and Owren generalised this concept and developed a framework for deriving linearly implicit energy-preserving multi-step methods for Hamiltonian PDEs with polynomial invariants [28]. A comparison of this approach and Kahan’s method applied to PDEs is given in [29]. Recently, more work has been put into developing linearly implicit energy-preserving schemes for Hamiltonian PDEs, e.g. the partitioned averaged vector field (PAVF) method [30] and schemes based on the invariant energy quadratization (IEQ) approach [31] or the multiple scalar auxiliary variables (MSAV) approach [32]. However, little attention has been given to linearly implicit local energy-preserving methods. To the best of the authors’ knowledge, the only existing method is one based on the IEQ approach, specific for the sine-Gordon equation [33]. In this paper, we use Kahan’s method to construct a linearly implicit method that preserves a discrete approximation to the local energy for multi-symplectic PDEs with a cubic energy function. This class is extensive, and includes among other PDEs the Korteweg–de Vries (KdV) equation, the Benjamin–Bona–Mahony (BBM) equation [34], Boussinesq-type systems [35], and the Camassa–Holm equation [36].

The rest of this paper is organized as follows. First, we give an overview of Kahan’s method and formulate it by using a polarised energy function. A brief introduction to multi-symplectic PDEs and their conservation laws are presented in Section 3. In Section 4, new linearly implicit local and global energy-preserving schemes are presented. Numerical examples for the KdV and Zakharov–Kuznetsov equations are given in Section 5, before we end the paper with some concluding remarks.

2 Kahan’s method

Consider an ODE system

[TABLE]

where $\hat{Q}(y)$ is an $\mathbb{R}^{M}$ valued quadratic form, $\hat{B}\in\mathbb{R}^{M\times M}$ is a symmetric constant matrix, and $\hat{c}\in\mathbb{R}^{M}$ is a constant vector. Kahan’s method is then given by

[TABLE]

where

[TABLE]

is the symmetric bilinear form obtained by polarisation of the quadratic form $\hat{Q}$ [24]. Polarisation, which maps a homogeneous polynomial function to a symmetric multi-linear form in more variables, was used to generalise Kahan’s method to higher degree polynomial vector fields in [37].

Suppose we restrict the problem (2.1) to be a Hamiltonian system on a Poisson vector space with a constant Poisson structure:

[TABLE]

where $A$ is a constant skew-symmetric matrix, and $H:\mathbb{R}^{M}\rightarrow\mathbb{R}$ is a cubic polynomial function. We first consider the Hamiltonian $H$ to be homogeneous. Then, following the result in Proposition $2.1$ of [37], Kahan’s method can be reformulated as

[TABLE]

where $\bar{H}(\cdot,\cdot,\cdot):\mathbb{R}^{M}\times\mathbb{R}^{M}\times\mathbb{R}^{M}\rightarrow\mathbb{R}$ is a symmetric $3$-tensor satisfying $\bar{H}(x,x,x)=H(x)$. Consider the $3$-tensor $\bar{H}(x,y,z)=x^{T}Q(y)z$, where $Q(y)=\frac{1}{6}\nabla^{2}H(y)$, with $\nabla^{2}H$ being the Hessian of $H$; then we can rewrite Kahan’s method (2.3) as

[TABLE]

where $\frac{\partial\bar{H}}{\partial x}$ denotes the partial derivative with respect to the first argument of $\bar{H}$.

Consider then the cases where the Hamiltonian in problem (2.2) is non-homogeneous, i.e. of the general form

[TABLE]

where $Q(y)$ is the linear part of $\nabla^{2}H(y)$ and thus a symmetric matrix whose elements are homogeneous linear polynomials, $B$ is the constant part of $\nabla^{2}H(y)$ and thus a symmetric constant matrix, $c$ is a constant vector and $d$ is a constant scalar. We follow the technique in [24], adding one variable to $y=(y_{1},\dots,y_{M})^{T}$ to get $\tilde{y}=(y_{0},y_{1},\dots,y_{M})^{T}$, extending $A$ to $\tilde{A}$ by adding a zero initial row and a zero initial column, considering a homogeneous function $\tilde{H}(\tilde{y})$ based on the non-homogeneous Hamiltonian $H(y)$ such that $\tilde{H}(\tilde{y})\lvert_{y_{0}=1}=H(y)$, and finally solving instead of (2.2) the equivalent, homogeneous cubic Hamiltonian problem

[TABLE]

with $y_{0}=1$. In this way we can still get the reformulation of Kahan’s method as (2.4) with

[TABLE]

Remark 1**.**

The $\mathbb{R}$-valued function $\bar{H}(x,y,z)$ in (2.6) has the following properties:

$\bar{H}(x,y,z)$* is symmetric111Denote the elements in $Q(y)$ by $q_{ij}y=\sum_{k}q_{ij}^{k}y_{k}$, where $q_{ij}^{k}$, $i,j,k=1,\cdots,M$, are scalars and $y_{k}$ is the $k$th element of $y$. We have that $q_{ij}^{k}$ satisfies $q_{ij}^{k}=q_{ki}^{j}=q_{jk}^{i}$ since $q_{ij}^{k}=\frac{1}{6}\frac{\partial^{3}\bar{H}}{\partial y_{i}\partial y_{j}\partial y_{k}}$, which is unchanged under any permutation of $i,j,k$. This provides the symmetry of $\bar{H}(x,y,z)$. w.r.t. $x$, $y$ and $z$,* 2. 2.

$\bar{H}(x,x,x)=H(x)$, 3. 3.

$\frac{\partial\bar{H}(x,y,z)}{\partial x}=Q(y)z+\frac{B(y+z)}{3}+\frac{c}{3}$* is symmetric w.r.t. $y$ and $z$.*

In this paper, we will use the form of Kahan’s method in (2.4) to prove the energy preservation of the proposed methods.

3 Conservation laws for multi-symplectic PDEs

Many PDEs, including all one-dimensional Hamiltonian PDEs, can be written on the multi-symplectic form

[TABLE]

where $K$, $L\in\mathbb{R}^{l\times l}$ are two constant skew-symmetric matrices and $S:\mathbb{R}^{l}\mapsto\mathbb{R}$ is a scalar-valued function. Following the results about multi-symplectic structure in [5], it can be shown that multi-symplectic PDEs satisfy the following local conservation laws [38]: the multi-symplectic conservation law

[TABLE]

the local energy conservation law (LECL)

[TABLE]

and the local momentum conservation law (LMCL)

[TABLE]

where $K_{+}$ and $L_{+}$ satisfy

[TABLE]

Decomposition of the matrices is done to make deduction of the conservation laws for energy and momentum more efficient [13, Section 12.3.1].

The multi-symplectic form (3.1) can also be generalised to problems in higher dimensional spaces. Consider $d$ spatial dimensions; based on the work by Bridges [5], a multi-symplectic PDE can then be written as

[TABLE]

where $K$, $L^{\alpha}\in\mathbb{R}^{l\times l}$ $(\alpha=1,\ldots,d)$ are constant skew-symmetric matrices and $S:\mathbb{R}^{l}\rightarrow\mathbb{R}$ is a smooth functional. Equation (3.3) has the following local energy conservation law:

[TABLE]

where $E(z)=S(z)+\sum_{\alpha=1}^{d}z_{\alpha}^{T}L^{\alpha}_{+}z$, $F^{\alpha}=-z_{t}^{T}L^{\alpha}_{+}z$, and $L^{\alpha}_{+}$ are splittings of $L^{\alpha}$ satisfying $L^{\alpha}=L^{\alpha}_{+}-({L^{\alpha}_{+}})^{T}$.

Say we have (3.3) defined on the spatial domain $\Omega\in\mathbb{R}^{d}$ with periodic boundary conditions. Integrating over the domain $\Omega$ on both sides of the equation (3.4) and using the periodic boundary condition then leads to the global energy conservation law for the multi-symplectic PDEs,

[TABLE]

where $\mathcal{E}(z)=\int_{\Omega}E(z)d\Omega$.

Example 1**.**

Korteweg–de Vries equation. Consider the KdV equation for modeling shallow water waves,

[TABLE]

where $\eta,\gamma\in\mathbb{R}$. Introducing the potential $\phi_{x}=u$, momenta $v=\gamma u_{x}$ and the variable $w=\gamma v_{x}\phi_{t}+\frac{\gamma^{2}u^{2}}{2}$ by the covariant Legendre transform from the Lagrangian, we obtain

[TABLE]

from which we find the multi-symplectic formulation (3.1) for the KdV equation with $z=(\phi,u,v,w)^{T}$, the Hamiltonian $S(z)=\frac{v^{2}}{2}-uw+\frac{\eta u^{3}}{6}$, and

[TABLE]

As for the conservation laws, there are many choices of $K_{+}$ and $L_{+}$, for example $K_{+}=\frac{K}{2},L_{+}=\frac{L}{2}$, or $K_{+}$ and $L_{+}$ being the upper triangular parts of $K$ and $L$, respectively.

Example 2**.**

Zakharov–Kuznetsov equation. Zakharov and Kuznetsov introduced in [39] a (2+1)-dimensional generalisation of the KdV equation which includes weak transverse variation,

[TABLE]

A multi-symplectification of this leads to a system (3.3) for two spatial dimensions,

[TABLE]

Following [11], we have that (3.8) is equivalent to a system of first-order PDEs,

[TABLE]

which is (3.9) with $z=(p,u,q,\phi,v,w)^{T}$, the Hamiltonian $S(z)=up-\frac{1}{2}(v^{2}+w^{2})-\frac{1}{6}u^{3}$, and the skew-symmetric matrices $K,L^{1},L^{2}$ whose only non-zero elements are

[TABLE]

4 New linearly implicit energy-preserving schemes

In [21], Gong, Cai and Wang present a scheme that preserves the local energy conservation law (3.4) of a one-dimensional multi-symplectic PDE, obtained by applying the midpoint rule in space and the averaged vector field (AVF) method in time. They also present schemes that preserve the global energy, but not (3.4), obtained by considering spatial discretizations that preserve the skew-symmetric property of the difference operator $\partial_{x}$. We build on their work by considering Kahan’s method for the discretization in time, ensuring linearly implicit schemes and also energy preservation.

To introduce our new schemes, we begin with some basic difference operators:

[TABLE]

The operators satisfy the following properties [18]:

All the operators commute with each other, e.g.

$\delta_{t}\delta_{x}v_{j}^{n}=\delta_{x}\delta_{t}v_{j}^{n},\quad\delta_{t}\mu_{x}v_{j}^{n}=\mu_{x}\delta_{t}v_{j}^{n},\quad\mu_{t}\delta_{x}v_{j}^{n}=\delta_{x}\mu_{t}v_{j}^{n}$. 2. 2.

They satisfy the discrete Leibniz rule

[TABLE]

Specifically,

[TABLE]

One can obtain a series of similar commutative equations and discrete Leibniz rules that are not presented here, but which are also crucial in the proofs of the preservation properties of the schemes to be introduced in the remainder of this section.

4.1 A local energy-preserving scheme for multi-symplectic PDEs

In this section, we apply the midpoint rule in space and Kahan’s method in time to construct a local energy-preserving method for multi-symplectic PDEs. Introducing the concept by first considering the one-dimensional system (3.1), we apply the midpoint rule in space to get

[TABLE]

Then applying Kahan’s method gives us the linearly implicit local energy-preserving (LILEP) scheme

[TABLE]

Here we consider $S$ of the form $S(y)=y^{T}Q(y)y+y^{T}By+c^{T}y+d$, as in (2.5), and accordingly $\bar{S}(x,y,z)$ of the form (2.6).

Theorem 1**.**

The scheme (4.1) satisfies the discrete local energy conservation law

[TABLE]

where

[TABLE]

Proof.

Taking the inner product with $\frac{1}{3}\delta_{t}\mu_{x}z_{j}^{n}$ on both sides of (4.1) and using the skew-symmetry of matrix $K$, we have

[TABLE]

Taking the inner product with $\frac{1}{3}\delta_{t}\mu_{x}z_{j}^{n+1}$ on both sides of (4.1), we get

[TABLE]

Taking the inner product with $\frac{1}{3}\delta_{t}\mu_{x}z_{j}^{n}$ on both sides of the scheme (4.1) for the next time step, we get

[TABLE]

Adding equations (4.4), (4.5) and (4.6) and using the skew-symmetry of matrix $K$, we obtain

[TABLE]

On the other hand, using the aforementioned commutative laws and discrete Leibniz rules for the operators, we can deduce

[TABLE]

Using the above relations (4.8), the fact that $L=L_{+}-L_{+}^{T}$ and the result (4.7), we obtain

[TABLE]

∎

Corollary 1**.**

For periodic boundary conditions $z(x+P,t)=z(x,t)$, the scheme (4.1) satisfies the discrete global energy conservation law

[TABLE]

where $\Delta x=P/M$ and $(\bar{E}_{L})_{j}^{n}$ is given by (4.3).

Proof.

With periodic boundary conditions, we get $\sum_{j=0}^{M-1}\delta_{x}(\bar{F}_{L})_{j}^{n}=0,$ and thus (4.9) follows from (4.2). ∎

The polarised global energy $\bar{\mathcal{E}}_{L}^{n}$ may be considered as a function of the solution in time step $n$ only, similarly to the modified Hamiltonian defined in Proposition 3 of [24].

Proposition 1**.**

With the solution $z^{n+1}$ found from $z^{n}$ by (4.1), the discrete global energy $\bar{\mathcal{E}}_{L}^{n}$ of (4.9) satisfies

[TABLE]

where

[TABLE]

while $z_{j}^{n+1}-z_{j}^{n}$ satisfies

[TABLE]

with $g_{L}(z_{j}^{n})=\nabla S(\mu_{x}z_{j}^{n})-L\delta_{x}z_{j}^{n}$ and $R_{L}(z_{j}^{n})=K\mu_{x}-\frac{\Delta t}{2}\nabla g_{L}(z_{j}^{n})$.

Proof.

Note that

[TABLE]

and

[TABLE]

Inserting this in (4.3), we get (4.10) from (4.9). Furthermore, observing that

[TABLE]

we may rewrite (4.1) as

[TABLE]

which is (4.12). ∎

Note that (4.11) is the discrete energy preserved by the fully implicit local energy-preserving method of [21]. Also, for methods based on the multi-symplectic structure, instead of solving for $z$ directly, the normal procedure is to eliminate the auxiliary variables from the scheme and get an equation for one variable $u$. Therefore we do not give an explicit expression for the modified energy in $z^{n}$. However, in Section 5, we present an explicit expression for the modified energy in $u^{n}$ when our scheme is applied to the KdV equation.

The results about the energy conservation for the LILEP method applied to one-dimensional multi-symplectic PDEs can be generalised to problems in spatial dimensions of any finite degree. Consider for example a $2$-dimensional multi-symplectic PDE

[TABLE]

for which we have the following corollary. This is presented without its proof, which is rather technical but similar to the proof of Theorem 1.

Corollary 2**.**

The scheme obtained by applying the midpoint rule in space and Kahan’s method in time to equation (4.13),

[TABLE]

where $j=0,\ldots,M_{x}-1$ and $k=0,\ldots,M_{y}-1$, satisfies the discrete local energy conservation law

[TABLE]

where

[TABLE]

4.2 Global energy-preserving methods for multi-symplectic PDEs

As shown in Section 3, Hamiltonian PDEs of the form (3.1) with periodic boundary conditions have global energy conservation which can be deduced from the local conservation law. On the other hand, the local conservation law is not inherent in the global conservation law. In this section, we will focus on giving a systematic method that preserves the global energy conservation law directly. We discretize $\partial_{x}$ with an antisymmetric differential matrix $D$ and get the semi-discretized variant of (3.1),

[TABLE]

where $z:=(z_{0},z_{1},\ldots,z_{M-1})^{T}\in\mathbb{R}^{M\times l}$ and $(Dz)_{j}=\sum_{k=0}^{M-1}D_{j,k}z_{k}$. We then apply Kahan’s method to (4.15) and obtain the linearly implicit global energy-preserving (LIGEP) scheme

[TABLE]

Define the polarised energy density by

[TABLE]

and we get the following result.

Theorem 2**.**

For periodic boundary conditions $z(x+P,t)=z(x,t)$, the scheme (4.16) satisfies the discrete global energy conservation law

[TABLE]

Proof.

Taking the inner product with $\frac{1}{3}\delta_{t}z_{j}^{n}$ on both sides of equation (4.16) and using the skew-symmetry of the matrix $K$, we get

[TABLE]

Taking the inner product with $\frac{1}{3}\delta_{t}z_{j}^{n+1}$ on both sides of (4.16), we get

[TABLE]

Furthermore, taking the inner product with $\frac{1}{3}\delta_{t}z_{j}^{n}$ on both sides of (4.16) for the next time step, we have

[TABLE]

Adding equations (4.19), (4.20) and (4.21), we get

[TABLE]

By using the commutative laws and discrete Leibniz rules,

[TABLE]

Based on the above equations (4.22) and (4.23), we obtain

[TABLE]

where

[TABLE]

Since $D$ is skew-symmetric and $G_{j,k}=G_{k,j}$, we get

[TABLE]

which implies that the discrete global energy conservation law $\bar{\mathcal{E}}^{n+1}=\bar{\mathcal{E}}^{n}$ is satisfied. ∎

The polarised energy $\bar{\mathcal{E}}$ preserved by (4.16) may also be expressed as a modification of the discrete energy

[TABLE]

which is preserved by the fully implicit global energy-preserving scheme of [21]. The proof of the following proposition is similar to the proof of Proposition 1, and hence omitted.

Proposition 2**.**

If the solution $z^{n+1}$ is found from $z^{n}$ by (4.16), the discrete global energy $\bar{\mathcal{E}}^{n}$ of (4.18) satisfies

[TABLE]

and $z_{j}^{n+1}-z_{j}^{n}$ satisfies

[TABLE]

where $g(z_{j}^{n})=\nabla S(z_{j}^{n})-L(Dz)_{j}^{n}$ and $R(z_{j}^{n})=K+\frac{\Delta t}{2}\nabla g(z_{j}^{n})$.

The above global conservation results can be generalised to multi-symplectic formulations in higher spatial dimensions, as demonstrated for the two-dimensional case by the following corollary, whose omitted proof is in the same vein as the proof of Theorem 2.

Corollary 3**.**

Discretizing $\partial_{x}$ and $\partial_{y}$ by skew-symmetric differential matrices $D_{x}$ and $D_{y}$ in equation (4.13) and then applying Kahan’s method to the semi-discrete system, one obtains the linearly implicit global energy-preserving (LIGEP) scheme

[TABLE]

where $j=0,\ldots,M_{x}-1$ and $k=0,\ldots,M_{y}-1$. For periodic boundary conditions $z(x+P_{x},y,t)=z(x,y,t)$, $z(x,y+P_{y},t)=z(x,y,t)$, the scheme (4.25) satisfies the discrete global energy conservation law

[TABLE]

where

[TABLE]

5 Numerical examples

In this section, we apply our proposed new linearly implicit energy-preserving schemes to the KdV equation and Zakharov–Kuznetsov equation, and compare them with fully implicit schemes. Among our reference methods are the methods introduced in [21], for which the local energy-preserving method is denoted by LEP, and the global energy-preserving method by GEP. These schemes are discretized in space the same way as our LILEP and LIGEP schemes, but the fully implicit AVF method is used for the time-stepping. For the GEP and LIGEP schemes, two different choices are considered for approximating the spatial derivative: the central difference operator $\delta^{c}_{x}$ defined by $\delta^{c}_{x}v_{j}^{n}:=\frac{1}{2}(\delta_{x}v_{j-1}^{n}+\delta_{x}v_{j}^{n})$ and the first order Fourier pseudospectral operator [11]. The latter results in the $M\times M$ matrix $D$, given explicitly by its elements

[TABLE]

evaluated on the domain $\left[0,P\right]$, where we assume $M$ even and periodic boundary conditions [40]. If $M$ is odd, we have instead

[TABLE]

The numerical results presented in this section are obtained from schemes implemented in MATLAB (2018b release), running on an early 2015 MacBook Pro with a dual-core 3.1 GHz Intel Core i7 processor and 16 GB of 1867 MHz DDR3 RAM. All fully implicit schemes are solved at each step by Newton’s method until $\lVert F(u^{n})\rVert_{2}<10^{-10}$. Linear systems are solved using the backslash operator of MATLAB. MATLAB and Python codes for the experiments are available at https://doi.org/10.5281/zenodo.3709463.

5.1 Korteweg–de Vries equation

Consider the multi-symplectic structure of the KdV equation as presented in Example 1. Applying the LILEP scheme (4.1) to (3.7), we obtain

[TABLE]

By eliminating the auxiliary varibles $\phi,v$ and $w$, we see that this is equivalent to

[TABLE]

Omitting the average operator $\mu_{t}$ gives us

[TABLE]

The polarised discrete energy preserved by this scheme is

[TABLE]

On the other hand, the discrete energy preserved by the LEP method of [21] is

[TABLE]

By Proposition 1 and elimination of the variables $\phi,v$ and $w$, (5.2) can be expressed as a modification of (5.3): we may rewrite (5.1) as

[TABLE]

where $(\mu_{x}u^{n})^{2}$ denotes the element-wise square of $\mu_{x}u^{n}$. Inserting this in (5.2), we get

[TABLE]

with

[TABLE]

where $\nabla\mathcal{E}_{L}^{n}$ means the gradient of $\mathcal{E}_{L}^{n}$ with respect to $u^{n}$, and $\zeta_{L}^{\prime}(u^{n})$ denotes the Jacobian matrix of $\zeta_{L}(u^{n})$.

Similarly for the LIGEP method (4.16); applying it to the the multi-symplectic KdV equations (3.7) and eliminating the auxiliary varibles $\phi,v$ and $w$, we obtain

[TABLE]

where $u^{n}u^{n+1}$ denotes element-wise multiplication of the vectors. Omitting the average operator $\mu_{t}$, we get

[TABLE]

The discrete global energy preserved by the GEP method is

[TABLE]

while the polarised discrete energy preserved by (5.4) is

[TABLE]

where $\zeta(u^{n})=-\gamma^{2}D^{3}u^{n}-\frac{\eta}{2}D(u^{n})^{2}$.

Test problem 1

In the first numerical experiment, we consider the problem introduced in [41] and then used by Zhao and Qin [42] and Ascher and McLachlan [43] to test various symplectic and multi-symplectic schemes: the KdV equation with $\gamma=0.022$, $\eta=1$, and initial value

[TABLE]

with $x\in[0,P]$, $P=2$. This problem is also considered in Example 3 of [21], where it is solved by implicit schemes that preserve local and/or global energy. As observed by Gong et al.  the global energy-preserving scheme (GEP) with the central difference operator used to approximate $\partial_{x}$ gives unsatisfactory results for this problem; we observed that the same is true for the LIGEP scheme. Therefore, the Fourier pseudospectral operator is used to approximate the spatial derivatives in the GEP and LIGEP schemes. This seems to result in more accurate solutions than the LEP and LILEP schemes for the same number of discretization points, but at a considerably higher computational cost, as seen from Table 1. Also from Table 1, we see an example of the advantage that can be gained by having a linearly implicit scheme instead of a fully implicit scheme. The different running times give an indication of the number of iterations necessary in Newton’s method to solve the fully implicit schemes for the particular cases. From Figure 1, we can conclude that our linearly implicit schemes give results close to their fully implicit counterparts introduced in [21], and that the different schemes converge to the same solution.

Compared to the schemes tested in [42, 43], our schemes do also perform well; see Figure 2, where we have plotted solutions by our schemes for the same discretization parameters used in Example 5.3 of [43]. The reference solution is found by the implicit midpoint scheme of [43] with $\theta=1$ and very fine discretization in space and time: $M=2000$ and $\Delta t=0.0001$. We observe that the LILEP scheme behaves similarly to the multi-symplectic box scheme of Arscher and McLachlan (see figures 3 and 4 in [43]), seemingly with the same superior stability for rough discretization in space and time. The LIGEP scheme, on the other hand, starts to blow up at around $t=1$ when $M=60$, $\Delta t=1/150$, but produces for $M=100$, $\Delta t=0.004$ a solution that is much closer to the correct solution than any of the schemes tested in [43] (see Figure 3 in that paper for comparison).

Test problem 2

To get quantitative results on the performance of our methods, we wish to study a problem with a known solution. For the KdV equation with $\gamma=1$, $\eta=6$, initial value $u_{0}(x)=\frac{1}{2}c\,\mathrm{sech}^{2}(-x+P/2)$ and periodic boundary conditions $u(x+P,t)=u(x,t)$, the exact solution is a soliton moving with a constant speed $c$ in the positive $x$-direction while keeping its initial shape. That is,

[TABLE]

In our numerical experiments, $c=4$ and $P=20$. For this problem, we have used the central difference operator to approximate $\partial_{x}$ in the GEP and LIGEP schemes, since it gives good results and yields considerably shorter computational time than if the pseudospectral operator is used. The proposed methods all show very good stability conditions when applied to this problem, as expected by methods conserving some invariant. The initial shape of the wave is well kept for long integration times, even when quite large step sizes in space and time are used; Figure 3 gives a good illustration of this. As in the previous example, we again observe that little is lost in accuracy by choosing linearly implicit over fully implicit time integration. A close inspection of Figure 3 also indicates that the local energy-preserving schemes preserve the shape of the wave better than the global energy-preserving schemes, while on the other hand, the GEP and LIGEP schemes are better than the LEP and LILEP schemes at preserving the speed of the wave. This is confirmed in Table 2 by measuring the shape error

[TABLE]

and phase error

[TABLE]

where $U^{N}$ is the numerical solution at end time $t$.

In Figure 4, we have plotted the computational time required to reach a certain accuracy in the global error for the different methods, both at time $t=0.5$ and at time $t=10$. We compare our methods to the fully implicit LEP and GEP schemes of [21], to a scheme based on discretizing the standard form (3.6) of the KdV equation in space and applying the AVF method in time (AVFM), and also to two of the schemes studied in [43]: the multi-symplectic box scheme (MSB) and the implicit midpoint scheme (IMP). Most notably we see from both plots in Figure 4 that the linearly implicit schemes perform better than the fully implicit schemes. Also, we see that at time $t=0.5$ the global error is lowest for the LILEP scheme, while at $t=10$ it is lowest for the LIGEP scheme. This is in accordance with the schemes’ phase and shape errors, which can be observed from Figure 3 and Table 2; with increasing time, the phase error becomes more dominant, and thus the scheme with the smallest phase error becomes increasingly advantageous.

Figure 5 illustrates how the different schemes preserve a discrete approximation to the energy to machine precision. That is, the linearly implicit schemes LILEP and LIGEP preserve exactly the discrete energies (5.2) and (5.6), respectively, while keeping the discrete energies (5.3) and (5.5), respectively, within some bound which depends on the discretization parameters. Likewise, the reverse is true for the fully implicit schemes. These observations fit well with our above results about the different discrete approximations to the energy: that for both the local energy preserving and the global energy preserving schemes, either discrete energy given can be seen as a modification of the other approximation. Finally, we have included plots in Figure 6 which confirm that our schemes are of second order in space and time.

5.2 Zakharov–Kuznetsov equation

Kahan’s method is previously shown to have nice properties when applied to integrable ODE systems [24, 25], and to perform well compared to other linearly implicit methods when applied to the KdV and Camassa–Holm equations [29], which are completely integrable PDEs. We wish to test our methods also on non-integrable systems, as well as on higher-dimensional problems. Therefore we consider the Zakharov–Kuznetsov equation, which is a non-integrable PDE [44, 45]. This two-dimensional generalisation of the KdV equation has a variety of applications, see e.g. [46] for a brief summary.

Applying the LILEP method (4.14) to the Zakharov–Kuznetsov equation (3.8) multi-symplectified as described in Example 2, we find

[TABLE]

Upon eliminating all variables except $u$, we are left with

[TABLE]

The operator $\mu_{t}$ is again superfluous. Hence we get the scheme

[TABLE]

This scheme preserves

[TABLE]

which is a two-step discrete approximation of the energy

[TABLE]

Similarly, applying the linearly implicit global energy-preserving method (4.25) to (3.10), we get the scheme

[TABLE]

which preserves the two-step discrete energy approximation

[TABLE]

Test problem

Taking a note from a numerical experiment performed in [11], we study the formation of cylindrical soliton pulses on the domain $\left[0,P\right]\times\left[0,P\right]$, $P=30$, following the initial condition

[TABLE]

where $\xi(y)$ is a random perturbation.

Upon trying the different schemes we can immediately conclude that the local energy-preserving schemes are superior for this problem when compared to the global energy-preserving schemes. The GEP and LIGEP schemes are too costly when the pseudospectral operator is used, and gives oscillatory behaviour in the $y$-direction when the central difference operator is used, unless the discretization in this direction is very fine. Although the global energy-preserving schemes with the central difference operator are slightly faster then the local energy-preserving schemes, as can be seen in Table 3, this is undermined by the cost of the extra discretization points needed to avoid oscillations in the former case. As was the case for the KdV problem, we see little difference between the linearly implicit schemes and their fully implicit counterparts. This can be seen in Figure 7, as can the oscillations in $y$-direction of the solution found by the GEP and LIGEP methods. The plots in Figure 7 can be compared to the plot in Figure 8, where the same problem is solved by the LILEP method using finer discretization in space and time. The initial random perturbation in $y$-direction over $75$ points is then transferred over to $225$ points using linear interpolation.

6 Concluding remarks

In this paper, we propose two types of linearly implicit methods with conservation properties for cubic invariants of multi-symplectic PDEs. The linearly implicit local energy-preserving (LILEP) method preserves a discrete approximation to the local energy conservation law, and by extension, the global energy whenever periodic boundary conditions are considered. The linearly implicit global energy-preserving (LIGEP) method preserves the global energy without inheriting the local preservation from the continuous system.

We test our methods on two PDEs: the one-dimensional, integrable Korteweg–de Vries (KdV) equation and the two-dimensional, non-integrable Zakharov–Kuznetsov equation. The numerical experiments confirm that the proposed methods are of second order both in space and time and that they preserve the expected local and global energy conservation laws. We have observed excellent stability properties for the LILEP scheme in particular, and very high accuracy in the LIGEP scheme even for quite coarse discretization when a Fourier pseudospectral operator is used to approximate the spatial derivative. Compared to the fully implicit methods of Gong et al. in [21], which was an inspiration for this paper, our methods show comparable wave profiles, global errors and energy errors, at a significantly lower computational cost. For two-dimensional problems, where fully implicit schemes quickly become very expensive to compute, the combination of local energy-preservation and a linearly implicit method seems to provide for a very competitive method.

Although we have only considered the preservation of cubic invariants in this paper, our schemes can be extended to preserve higher order polynomials by the polarisation techniques for generalising Kahan’s method suggested in [37]. This would result in $(p-2)$-step methods for preservation of a discrete $p$-order polynomial invariant. Although the idea behind this is clear, obtaining the concrete schemes is not straightforward, and thus we leave it for future research.

Acknowledgements

This work was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 691070. The authors wish to express gratitude to Elena Celledoni and Brynjulf Owren for constructive discussions and helpful suggestions during our work on this paper, and to Benjamin Tapley for helping with the language.