A New Technique for Preserving Conservation Laws

Consider a system of q partial differential equations (PDEs),

where \({\mathcal {A}}\) is a row vector, \({\mathbf {u}}\) has components \(u^\alpha ,\ \alpha =1,\dots ,q\), and square brackets around a differentiable expression denote the expression and finitely many of its derivatives.^{Footnote 1} We assume that (1) is totally nondegenerate (see [18]). A local conservation law is a divergence expression,

that vanishes on all solutions of (1). The functions F and G are the flux and the density of the conservation law, respectively, and \(D_x\) and \(D_t\) denote the total derivatives with respect to x and t, respectively. The conservation law (2) is in characteristic form if there exists a column vector \({\mathcal {Q}}\) such that

in which case \({\mathcal {Q}}\) is called the characteristic. The space of total divergences forms the kernel of the Euler operator, \({\mathcal {E}}\), whose \(\alpha \)th entry is

Hence

if and only if there exists \({\mathbf {F}}\) such that (3) holds. These results generalize immediately to PDE systems with more than two independent variables.

The literature on the numerical solution of PDEs is rich in numerical methods that preserve global invariants, but there are relatively few results on the preservation of local conservation laws. Arguably, local conservation laws are more necessary: they hold throughout the domain, apply to the set of all solutions, and provide much stronger constraints than are needed to preserve the corresponding global invariants. Moreover, when the domain and boundary conditions are suitable, conservation of such invariants is automatically achieved.

A new approach for developing finite difference schemes that preserve conservation laws of (1) was introduced recently in [9]. It exploits the fact that discrete conservation laws form the kernel of a discrete version of the Euler operator (4). Discretizations of the PDE (1) having discrete versions of the desired conservation laws are obtained by requiring that a discrete version of condition (5) is satisfied. This requires the symbolic solution of a large system of nonlinear equations that is impractical in general. The complexity of the symbolic calculations can be reduced by introducing compactness requirements on the schemes, and this direct approach has been applied to a range of PDEs with different structure in [9,10,11]. However, the direct approach is greatly limited by the capacity of symbolic computation; it has been applied only to second-order approximations of PDEs with two independent variables.

In this paper, we modify the approach in [9] by finding semidiscretizations of (1) that preserve semidiscrete local conservation laws. The reduction to one discrete space dimension significantly reduces the complexity of the computations, to the point that the determining system can be solved easily without introducing any restrictions on the form of the schemes. After this, a suitable integrator in time needs to be chosen to create a fully discrete scheme; this depends on the form of the conservation laws that one aims to preserve.

If the PDE is equipped with conservative boundary conditions, it is known that the quadratic invariants of its space discretizations are preserved by symplectic Runge–Kutta methods [4, 8, 20]. In this paper, we extend this result to prove that if G is quadratic in [u], then any symplectic Runge–Kutta method preserves the conservation law (2) locally, regardless of the boundary conditions.

There are various results on local conservation for Hamiltonian PDEs,

where \([{\mathbf {u}}]_x\) denotes \({\mathbf {u}}\) and its spatial derivatives only, \({\mathcal {D}}\) is a skew-adjoint operator that satisfies the Jacobi identity, and H is the Hamiltonian function.

Multisymplectic schemes [3] and their generalizations [21] can preserve local conservation laws with quadratic flux and density. Requiring the flux to be quadratic is, however, a strong constraint that is not satisfied by local momentum conservation laws of many important equations in physics such as the nonlinear Schrödinger (NLS) equation, the Korteweg–de Vries (KdV) equation, the Benjamin–Bona–Mahony (BBM) equation, the modified Korteweg–de Vries (mKdV) equation, and the Boussinesq equation. The strategy introduced in this paper does not suffer from this restriction, as no assumption is needed about the flux, so it can be applied to the preservation of these conservation laws as well.

Another popular approach is to use a discrete gradient method for the time integration of (6). These are obtained from a semidiscretization of H and a skew-adjoint discretization of \({\mathcal {D}}\), and preserve a discrete conservation law of the energy [17]. One widely used discrete gradient method is the average vector field (AVF) method [5, 6, 19]. We show that the AVF method yields the local conservation law of the Hamiltonian under constraints on the discretization of \({\mathcal {D}}\) that are milder than skew-adjointness. Consequently, conservation of the local Hamiltonian can be achieved for a larger class of discretizations.

Although the discussion so far has focused on PDEs with two independent variables, the approach of discretizing one dimension at a time works equally well for PDEs on higher-dimensional spaces. We discuss this and give an illustration.

The paper is organized as follows. Section 2 introduces a method for obtaining conservative spatial semidiscretizations. Section 3 focuses on time integration. In particular, we show the following.

• Conservation laws with quadratic density (without any assumption on the flux) are preserved by any symplectic method in time locally and independently of the boundary conditions. Conservation laws for mass, charge and momentum are typically in this class.

• For Hamiltonian PDEs, the AVF method preserves the local semidiscrete conservation law of the energy for a wide class of semidiscretizations that generalizes the result in [17].

• For other types of conservation law, fully discrete methods can be found by introducing relatively few parameters and fixing them by requiring that the conservation law is in the kernel of a fully discrete Euler operator. This approach can be iterated for dimensions, by using a sequence of semidiscrete and discrete Euler operators.

In Sect. 4, we apply this new approach to the Boussinesq equation and introduce methods of order two and four that preserve three conservation laws. Section 5 describes numerical benchmark tests, including evidence of stability and comparison with other methods from the literature. In Sect. 6, we apply the new technique to a two-dimensional PDE, the potential Kadomtsev–Petviashvili (pKP) equation and introduce two families of schemes that preserve two conservation laws. Finally, we draw some conclusions in Sect. 7.

We begin with a regular spatial grid. The stencil consists of \(M=B-A+1\) nodes,

where \(x_0\) is a generic grid point; let \({\mathbf {x}}\) denote the vector of the nodes. The forward shift operator \(S_m\) acts as follows on any semidiscrete function f:

the forward difference, forward average, and centred difference operators are

respectively, where I is the identity operator. The semidiscretizations of \(u^\alpha \!(x,\!t)\) are given by the column vector \({\mathbf {U}}\in {\mathbb {R}}^{Mq}\) with the \((m+\alpha M-B)\)th entry

The semidiscrete problem is

where here and henceforth tildes represent approximations to the corresponding continuous quantities, and square brackets around a semidiscrete expression denote the expression and a finite number of its time derivatives.

A semidiscrete conservation law of (9) is a semidiscrete divergence,

such that

The functions \({\widetilde{F}}\) and \({\widetilde{G}}\) are the semidiscrete flux and density of the conservation law (10), respectively. Similarly, as in the continuous case, we say that (10) is in characteristic form if there exists \({\widetilde{{\mathcal {Q}}}}={\widetilde{{\mathcal {Q}}}}({\mathbf {x}},t,[{\mathbf {U}}])\), called the characteristic, such that

The following result is crucial for obtaining semidiscretizations that preserve conservation laws (see [14, 18] for analogous results in the continuous and totally discrete setting, respectively).

Theorem 1

The kernel of the semidiscrete Euler operator \({\mathsf {E}}_{\mathbf {U}}\), whose \(\alpha \)th entry is

is the space of semidiscrete divergences (10).

Proof

Let \(L=L({\mathbf {x}},t,[{\mathbf {U}}])\) such that \({\mathsf {E}}_{\mathbf {U}}(L)={\mathbf {0}}\), and consider the derivative

Integrating by parts yields

for some functions \({\hat{F}}={\hat{F}}({\mathbf {x}},t,\varepsilon ,[{\mathbf {U}}])\) and \({\hat{G}}={\hat{G}}({\mathbf {x}},t,\varepsilon ,[{\mathbf {U}}])\) whose precise expression is not of importance. Substituting this into (11) and integrating over \(\varepsilon \in [0,1]\) shows that L is a semidiscrete divergence.

If L is of the form (10), \({\mathsf {E}}_{\mathbf {U}}(L)={\mathbf {0}}\) follows from the linearity of the Euler operator and from the fact that for any k (see e.g. [12]),

\(\square \)

Based on the result in Theorem 1, the approach used in [9,10,11] to preserve fully discrete conservation laws, is adapted here to the preservation of semidiscrete conservation laws of (1) with characteristics \({\mathcal {Q}}_\ell \), as follows:

1. 1.

   Select a stencil that is large enough to support generic semidiscretizations for \({\mathcal {A}}\) and every \({\mathcal {Q}}_\ell \), having the desired order of accuracy. These approximations depend on a number of free parameters to be determined.

2. 2.

   Find some of the parameters by imposing consistency, up to the desired order of accuracy, p.

3. 3.

   Use symbolic algebra to determine the values of the free parameters that satisfy

   $$\begin{aligned} {\mathsf {E}}_{\mathbf {U}}({\widetilde{{\mathcal {A}}}}{\widetilde{{\mathcal {Q}}}}_\ell )={\mathbf {0}}, \end{aligned}$$

   (12)

   for \(\ell =1\). This guarantees that the first conservation law is locally preserved. As both \({\widetilde{{\mathcal {A}}}}\) and \({\widetilde{{\mathcal {Q}}}}_\ell \) are accurate to order p, the discrete conservation law has the same order of accuracy.

4. 4.

   Iterate the previous step, replacing \({\widetilde{{\mathcal {Q}}}}_1\) with \({\widetilde{{\mathcal {Q}}}}_\ell \), to obtain further constraints on the parameters. If (12) has no solution for some \(\ell \), the corresponding conservation law cannot be preserved without violating one of the previous conservation laws. Typically, the more complicated a conservation law is, the more parameters need to be fixed to preserve it.

Remark 1

It might seem appealing to identify a set of conservation laws that one wishes to preserve and use brute force symbolic computation to solve all constraints simultaneously. (This was our approach initially.) However, this takes far longer than the sequential approach and commonly comes up with a null result, with no indication as to which subsets of conservation laws can be preserved. The sequential algorithm above enables the user to decide which conservation laws should be prioritized. At each iteration, the computation is simplified by the fact that some parameters have already been fixed.

Remark 2

If the algorithm above does not produce any schemes for a given stencil, one could try preserving the same conservation laws using a wider stencil. However, the wider the stencil is, the more the computational cost increases. Moreover, if one finds a conservative semidiscretization, a time integrator that preserves all the conservation laws is also needed. For example, in the next section we prove that some known time integrators preserve conservation laws whose density is either quadratic, or is a Hamiltonian, but to the best of our knowledge there are no methods that preserve both of these types.

We begin by considering one-step time integrators. For fully discrete schemes, the stencil is

and the forward shift operators in space and time are

respectively. The forward difference and forward average operators in space are defined by (8), and the corresponding operators in time are

Let \({\mathbf {u}}_{n}\in {\mathbb {R}}^{Mq}\) be the column vector whose \((m+\alpha M-B)\)th entry is

and let \({\mathbf {u}}_{m,n}\in {\mathbb {R}}^{q}\) be the column vector with entries

3.1 Conservation Laws with Quadratic Density

Here, attention is restricted to PDEs of the form

where \({\mathbf {g}}\) is linear homogeneous in \([{\mathbf {u}}]_x\); these include Hamiltonian PDEs. Consider a conservation law of (13) of the form

where the density, \(G_2\), is a polynomial of degree two in \([u]_x\). (Without loss of generality, assume that no terms in \(G_2\) depend on x only.) For many differential problems of importance in physics (such as KdV, NLS and BBM equations), the conservation laws of mass (or charge) and momentum are of the form (14) with linear and quadratic density, respectively.

Let \(P({\mathbf {x}})\) be an invertible operator such that

and \({\widetilde{{\mathbf {h}}}}({\mathbf {x}},t,{\mathbf {U}})\) are two column vectors whose \((m+\alpha M-B)\)th entry is a spatial discretization of the \(\alpha \)th component of \({\mathbf {g}}(x,[{\mathbf {u}}]_x)\) and \({\mathbf {h}}(x,t,[{\mathbf {u}}]_x)\) at \(x_m\), respectively. Let

be a semidiscretization of the PDE (13) with the following approximation to the conservation law (14):

Such a semidiscretization can be obtained using the technique in Sect. 2. The flux and density of (16) have the general form

where \(S({\mathbf {x}})=S({\mathbf {x}})^T\in {\mathbb {R}}^{Mq\times Mq}\) and \({\mathbf {w}}({\mathbf {x}})\in {\mathbb {R}}^{Mq}\) is a column vector.

The following theorem shows that symplectic Runge–Kutta methods preserve local conservation laws with quadratic density. The proof adapts Calvo, Iserles and Zanna’s proof that such methods preserve quadratic invariants of systems of ODEs [4], to take contributions from the flux into account.

Theorem 2

The solution of any symplectic Runge–Kutta method applied to (15) satisfies a discrete version of (16).

Proof

The conservation law (16) amounts to

Solving (15) using a s-stage symplectic Runge–Kutta method,

with internal stages

we obtain \({\mathbf {u}}_n=P(x){{\widetilde{{\mathbf {g}}}}}_n\). Moreover,

Using (19) to eliminate \({{\widetilde{{\mathbf {g}}}}}_n\) from the first sum and rearranging, gives

The condition of symplecticity,

and (17) give

which is an approximation of (14). \(\square \)

Remark 3

Multisymplectic methods preserve conservation laws whose density and flux are both quadratic. By contrast, Theorem 2 applies to all conservation laws that have quadratic density. As no assumption is needed on the flux, a larger class of conservation laws can be preserved.

The following results follow directly from the proof of Theorem 2.

Corollary 1

Any Runge–Kutta method preserves semidiscrete local conservation laws whose density is linear in \([{\mathbf {u}}]_x\).

Corollary 2

The symplectic implicit midpoint method (defined by (18)–(19) with \(s=1, b_1=1,\) and \(a_{1,1}=c_1=1/2\)) applied to (15) preserves the conservation law

3.2 Conservation Law for the Hamiltonian

We consider here the system of Hamiltonian PDEs (6) defined by a Hamiltonian function H on a domain with periodic boundary conditions. This assumption is introduced only for simplicity: the preservation of conservation laws is local and therefore independent of the specific boundary conditions assigned to the differential problem. The following local conservation law for the energy is satisfied by all solutions of (6):

with

Among the best-known energy-conserving discrete gradient methods is the average vector field (AVF) method [19]. We can use this in two different ways, depending on the number of points, M, in the stencil in (7). If M is odd, let \(A=-B\) so that the stencil is centred on \(x_0\); we denote a semidiscretization of \(H([{\mathbf {u}}]_x)\) on such a stencil by \({\widetilde{H}}({\mathbf {U}})\). The AVF method approximates (6) by

If M is even, let \(A=1-B\) so that the stencil is centred at the midpoint of \(x_0\) and \(x_1\). Denote a semidiscretization of \(H([{\mathbf {u}}]_x)\) as \({\widetilde{H}}(\mu _m{\mathbf {U}})\). We define the AVF method on such a stencil to be

McLachlan and Quispel proved in [17] that discrete gradient methods preserve a discrete version of (20) provided that \(\widetilde{{\mathcal {D}}}\) is a skew-adjoint approximation of \({\mathcal {D}}\). The following theorem proves that the AVF method (21) preserves the local conservation law for the energy (20) under a milder assumption.

Theorem 3

The AVF methods (21) and (22) preserve a discrete energy conservation law if there exists a function f defined on the stencil such that

or

respectively.

Proof

Equation (23) yields a discrete energy conservation law for (21), namely

where

Similarly, from (24), the conservation law preserved by (22) is

where

\(\square \) \(\square \)

Remark 4

Theorem 3 holds true in particular when \(\widetilde{{\mathcal {D}}}\) is skew-adjoint.

Remark 5

Condition (23) is satisfied if and only if

Similarly, condition (24) holds true if and only if

Conditions (25) and (26) provide simple practical tests for analysing whether the assumptions of Theorem 3 are satisfied by a given scheme.

Remark 6

The following Hamiltonian-preserving schemes can be obtained from (21) or (22), where the operator \(\widetilde{\mathcal D}\) is not skew-adjoint but satisfies (23) or (24), respectively:

• EC\(_8\) and the family of schemes MC\(_8\) for KdV in [9],

• EC\(_8\)(0) (which preserves a quartic density) and MC\(_8\)(0) for mKdV in [10],

• EC\(_6\) for BBM in [11].

3.3 Conservation Laws of Other Types and Multidimensional Domains

Fully discrete methods that preserve other types of conservation law can be obtained again by using an Euler operator. For simplicity, we restrict the discussion to the case of second-order schemes for a PDE with polynomial nonlinearity. These can be obtained by following the steps below:

1. 1.

   Let \(P=\prod _{i=1}^r L_i\) be a polynomial of degree r in the semidiscretization, where without loss of generality \(L_i\) is a linear approximation of a single monomial factor in the continuous counterpart of P. Therefore, \(L_i\) can depend on either \({\mathbf {U}}\) or \({\mathbf {U}}_t\). Assuming that the stencil has N points in the time dimension, discretize P using

   $$\begin{aligned} \sum _{j=0}^{N^r-1} \prod _{i=1}^r \alpha _j L_i({\mathbf {u}}_{i_j}), \qquad \alpha _j\in {\mathbb {R}}, \end{aligned}$$

   (27)

   where \(i_j\) is the ith digit of the representation of \(N^r-1-j\) in base N, ordered from right to left, and setting \(i_j=0\) if \(N^r-1-j\) has less than i digits. By varying \(i_j\) one obtains all possible combinations of N digits of length r. At this stage, we assume that the coefficients \(\alpha _j\) only satisfy the requirements for consistency. For example, in a one-step method linear quantities are approximated by

   $$\begin{aligned} L_1({\mathbf {U}})=L_1(\mu _n{\mathbf {u}}_0),\qquad L_1({\mathbf {U}}_t)=L_1(D_n{\mathbf {u}}_0), \end{aligned}$$

   and quadratic quantities by

   $$\begin{aligned}L_1L_2\approx&\, \alpha _0L_1({\mathbf {u}}_1)L_2({\mathbf {u}}_1)+\alpha _1L_1({\mathbf {u}}_0)L_2({\mathbf {u}}_1)\\&+\alpha _2L_1({\mathbf {u}}_1)L_2({\mathbf {u}}_0)+\alpha _3L_1({\mathbf {u}}_0)L_2({\mathbf {u}}_0). \end{aligned}$$
2. 2.

   The values of the parameters \(\alpha _j\) are obtained by solving

   $$\begin{aligned} {\mathsf {E}}_{{\mathbf {u}}}({\widetilde{{\mathcal {A}}}}{\widetilde{{\mathcal {Q}}}}_\ell )={\mathbf {0}}, \qquad \ell =1,2,\ldots , \end{aligned}$$

   (28)

   where \({\widetilde{{\mathcal {A}}}}\) and \({\widetilde{{\mathcal {Q}}}}_\ell \) are the approximations of the PDE and the characteristic obtained after steps 1 and 2, and \({\mathsf {E}}_{{\mathbf {u}}}\) is the difference Euler operator,

   $$\begin{aligned} {\mathsf {E}}_{{\mathbf {u}}}=\sum _{i,j}S_m^{-i}S_n^{-j}\frac{\partial }{\partial u_{i,j}}, \end{aligned}$$

   (29)

   whose kernel consists of all fully discrete conservation laws [14].

Remark 7

In total, net of consistency requirements, for each monomial of degree r one needs:

• \(\left( \begin{array}{c} M+r-1\\ r\end{array}\right) \) parameters for the semidiscretization. After solving (12), only a few of these will still be undetermined.

• \(N^r\) new parameters for the full discretization (27), to be determined by solving (28).

By contrast, if one searches directly for all full discretizations without semidiscretizing first, the strategy in [9] introduces \(\left( \begin{array}{c} NM+r-1\\ r\end{array}\right) \) variables for each monomial of degree r. This in general yields a huge nonlinear system whose solution is impractical.

The algorithm above can be iterated, to preserve conservation laws for PDEs with more than two independent variables, by discretizing a single variable at each iteration. At the kth iteration the Euler operator in (29) is replaced by

whose kernel consists in the space of conservation laws of the form

The proof is similar to that of Theorem 1.

We consider here the (Good) Boussinesq equation

cast as a system of two PDEs

This system can be written in Hamiltonian form,

with

System (31) has infinitely many independent conservation laws [1]. The first four are \(D_xF_i+D_tG_i=0\), with

with characteristics

respectively. When the boundary conditions are conservative (e.g. periodic), integrating in space (32)–(35) gives the preservation of the following invariants,

here \(I_3\) and \(I_4\) are the global momentum and the global energy, respectively.

4.1 Conservative Methods for the Boussinesq Equation

We look for semidiscretizations of (31) of the form

Solutions of (38) satisfy semidiscrete versions of the conservation laws (32) and (33) with \({\widetilde{{\mathcal {Q}}}}_1={\mathcal {Q}}_1\) and \({\widetilde{{\mathcal {Q}}}}_2={\mathcal {Q}}_2\). The linear and quadratic terms in (31) and (36) are approximated as

respectively, where \(Z_i\in \{U_i, V_i, D_t U_i, D_t V_i\}\), and the coefficients \(\alpha _i\) and \(\beta _{i,j}\) are chosen by requiring the desired order of accuracy and the preservation of the conservation law of either the momentum (34) or the energy (35), according to the strategy outlined in Sect. 2. We have not found any semidiscrete scheme that preserves both of these conservation laws, as the constraints on the approximations of the nonlinear terms that we have obtained from (12) are not compatible with each other.

In the following, we present the components of (38) and the characteristic \({\widetilde{{\mathcal {Q}}}}_{3/4}\) and density \({\widetilde{G}}_{3/4}\) of the remaining preserved conservation law. The corresponding flux \({\widetilde{F}}_{3/4}\) does not contribute to our global error estimates; in most cases, it has many terms and gives little insight, so we omit this.

Fully discrete schemes that preserve the local momentum or the local energy are then obtained by applying Gauss–Legendre method or AVF method, respectively. As these are Runge–Kutta methods, the conservation laws (32) and (33) are preserved as a consequence of Corollary 1.

Second-order schemes

Here, we introduce new families of second-order schemes that preserve three conservation laws. The stencil in space consists of four points, and we set \(A=-1\) and \(B=2\) in (7). All the free parameters in the formulae below \((\alpha ,\beta ,\gamma ,\xi )\) are \(O(\varDelta x^2).\) Free parameters corresponding to higher-degree perturbations are set equal to zero as their contribution is negligible.

Momentum-conserving schemes

We have obtained six families of semidiscretizations that preserve the conservation law (34), split by two different forms of the characteristic and the three parameter values \(s\in \{0,1/3,1/2\}\).

The first three families are

Three remaining are given by \({\widetilde{F}}_1, {\widetilde{F}}_2\) and \({\widetilde{G}}_2\) as in (39) together with

In the numerical tests section, we limit our investigations to the semidiscretizations obtained from (38) with (39) and \(s=\alpha =\beta =\xi =0\), \(\gamma =\lambda _1\varDelta x^2\); we use MC\(_2(\lambda _1)\) to denote the family of finite difference schemes obtained by using the symplectic implicit midpoint method (Gauss–Legendre method of order two) to discretize in time. The iterative technique described in Sect. 3.3 also finds these schemes, but no others.

Energy conserving schemes

There is only one family of semidiscretizations of the form (38) that preserves the local conservation law of the energy. For this family, \({\widetilde{F}}_1\) and \({\widetilde{G}}_2\) are defined as in (39) and

The resulting system of ODEs can be written in the form

with

The operator \(\widetilde{{\mathcal {D}}}\) is not skew-adjoint, but satisfies (24). Therefore, by applying the AVF method (22) to (41) we obtain a family of fully discrete schemes that preserve the local conservation law of the energy. We use EC\(_2(\lambda _2)\) to denote the schemes with \(\alpha =\gamma =0\) and \(\beta =\lambda _2\varDelta x^2\). Again, the iterative technique from Sect. 3.3 yields only these schemes.

Fourth-order schemes

The families of fourth-order schemes introduced here preserve three conservation laws and depend on free parameters (\(\alpha ,\beta ,\gamma ,\xi \)) that are all \(O(\varDelta x^4)\). Parameters introducing only perturbations of higher degree are set equal to zero.

For each of the semidiscretizations introduced in this section, the discrete fluxes \({\widetilde{F}}_j\) are second-order accurate, but \(D_m {{\widetilde{F}}}_j\) and the three preserved conservation laws,

are approximated with fourth-order accuracy.

Momentum-conserving schemes

On a spatial stencil with six points (\(A=-2\) and \(B=3\) in (7)), there are two families of semidiscretizations that preserve the local momentum conservation law. Let

The first family of semidiscretizations and their preserved conservation laws is given by

The second family has \({\widetilde{F}}_1,{\widetilde{F}}_2\) and \({\widetilde{G}}_2\) as in (42), together with

We use MC\(_4(\lambda _3)\) to denote the schemes obtained by applying the Gauss–Legendre method of order four to (42), with \(\alpha =\beta =\gamma =0\) and \(\xi =\lambda _3\varDelta x^4\).

Energy-conserving schemes

On the most compact six-point stencil, there are no semidiscretizations of the form (38) that preserve the local conservation law for energy. However, a seven-point stencil (\(B=-A=3\) in (7)) is more fruitful. For \(n\in {\mathbb {N}}\), let

and define the operators

and the functions

The family of semidiscretizations and conservation laws is as follows:

The systems of ODEs defined by (38) with (43) amount to

with

The operator \(\widetilde{{\mathcal {D}}}\), although not skew-adjoint, satisfies (23). We use EC\(_4(\lambda _4)\) to denote the family of schemes obtained by applying the AVF method of order four (see [19]) to (44) with \(\alpha =\lambda _4\varDelta x^4\) and \(\beta =\gamma =0\).

In this section, we solve a couple of benchmark problems to show the effectiveness and conservation properties of the numerical methods in Sect. 4.

The results are compared with the following second-order structure-preserving methods:

• The multisymplectic scheme for (30),

  $$\begin{aligned} \text {PS}\equiv&\,D_n^2\mu _m^4u_{-2,-1}-D_m^2\mu _m^2\mu _n^2u_{-2,-1}-D_m^2\mu _m\mu _n(\mu _m\mu _n u_{-2,-1})^2\\&+D_m^4\mu _n^2u_{-2,-1}=0, \end{aligned}$$

  developed in [22] and equivalent to the well-known Preissmann scheme.

• The symplectic scheme for (31) in [7],

  $$\begin{aligned} \text {MP}\equiv \big (&\,D_n u_{0,0}-D^{(c)}_m\mu _n v_{0,0},\\ \nonumber&D_n v_{0,0}+D^{(c)}_m(D_m^2\mu _n u_{-1,0}- (\mu _n u_{0,0})^2-\mu _n u_{0,0} \,\big )={\mathbf {0}}, \end{aligned}$$

  obtained by applying the midpoint rule to a suitable spatial discretization.

• The energy-conserving scheme for (31) in [16],

  $$\begin{aligned} \text {DVD}\equiv \big (&\,D_n u_{0,0}-D_m\mu _n v_{0,0},\, D_n v_{0,0}\\ {}&\!+\!D_m(D_m^2\mu _n u_{-2,0}\!-\!\tfrac{1}{3} (u_{-1,0}^2\!+\!u_{-1,0}u_{-1,1}\!+\!u_{-1,1}^2) \!-\!\mu _n u_{-1,0} )\big )\!=\!{\mathbf {0}}, \end{aligned}$$

  obtained using a discrete variational derivative method. This scheme can be obtained also by applying the AVF method to the Hamiltonian system of ODEs defined by

  $$\begin{aligned} \widetilde{{\mathcal {D}}}=\left( \begin{array}{cc} 0 &{} \quad D_mS_m^{-1}\\ D_m &{} \quad 0\end{array}\right) , \qquad H=\tfrac{1}{2}(U_{0}^2+V_{0}^2+\mu _m((D_mU_{-1})^2)+\tfrac{1}{3}U_{0}^3. \end{aligned}$$

To the best of our knowledge, there are no schemes in the literature for the Boussinesq equation that are fourth-order accurate in both space and time. Therefore, we compare the performance of the fourth-order schemes in Sect. 4 with the following finite difference scheme for (30) introduced in [13]:

The scheme FD\(_4\) is fourth-order accurate in space and second-order accurate in time, so to have a fair comparison we will use this scheme with a time step equal to \(\varDelta t^2\).

We consider \((x,t)\in \varOmega \equiv [a,b]\times [0,T]\) and periodic boundary conditions. We introduce on \(\varOmega \) a grid with \(I+1\) nodes, \(x_i\), in space and \(J+1\) nodes, \(t_j\), in time. Henceforth subscripts denote shifts from the point \((x_0,t_0)=(a,0)\) (e.g. \(u_{i,j}\simeq u(a+i\varDelta x,j\varDelta t)\)).

As the computational time is similar for all the schemes of the same order of accuracy, our comparisons are based on the error in the solution at the final time \(t = T\), evaluated as

where \(\Vert \cdot \Vert \) denotes the Euclidean norm. We also compare the errors in the global invariants (37) defined as

where \(\alpha =1,2,3,4.\) For the methods introduced in this paper, \({\widetilde{G}}_\alpha \) is given in Sect. 4. For all the other schemes, we set

Single soliton

For the first problem, we set \(\varOmega =[-60,60]\times [0,25]\) and the initial conditions given by the single soliton solution over \({\mathbb {R}}\),

where \(c=\sqrt{1-p^2}.\) We choose

Single soliton problem: plot of \(\pi (\varDelta x)\) showing convergence of the schemes

Full size image

We first examine the convergence and stability of the schemes found in Sect. 4, setting all parameters to zero. The order of convergence at various step sizes is measured by

and error\(_k\) denotes the error obtained from (45) with \(\varDelta x=\varDelta t=k/10\), for \(k=1,2,\ldots ,9\). The results in Table 1 and Fig. 1 show that all the methods tend to the exact solution with the expected maximum order of convergence as the grid is refined and are stable also for the largest step sizes.

Table 2 shows the error in the conservation laws and the solution for the different methods with \(\varDelta x=\varDelta t=0.5\). For this problem, the values of the free parameters that minimize the error in the solution of MC\(_2(\lambda _1)\), EC\(_2(\lambda _2)\), MC\(_4(\lambda _3)\), EC\(_4(\lambda _4)\) are \(\lambda _1=-0.21,\) \(\lambda _2=-0.20\), \(\lambda _3=0.06\), \(\lambda _4=0.03\). Such optimization is easy given that the solution is known, but is not currently feasible more generally. Therefore, the results obtained by setting the above parameters to zero are shown for comparison. This benchmark test illustrates that:

• All schemes preserve the first conservation law, but only those based on the formulation (31) preserve the second conservation law;

• The schemes introduced in this paper preserve three conservation laws up to machine accuracy;

• The new second-order schemes compare favourably with the methods from the literature in terms of accuracy in both the solution and the invariants that are not exactly conserved, even without optimizing the parameters;

• Choosing the optimal value of the free parameter, the error in the solution is roughly four times smaller (or more) than any other second-order method;

• The fourth-order methods are all more much accurate than FD\(_4\) with time step \(\varDelta t=0.5^2\), even for nonoptimal parameters.

Top: initial condition (dashed line) and solution of MC\(_2(-0.21)\) at \(T=25\) (solid line). Bottom: comparison of different schemes around the top of the soliton

Full size image

Single soliton: error of fourth-order schemes

Full size image

In Fig. 2, the upper plot shows the solution obtained by the most accurate scheme, MC\(_2(-0.21)\). The motion of the soliton does not produce any spurious oscillations. These can be seen (with amplitude of about \(10^{-2}\)) in the solutions of MP and DVD.

The lower plot shows the exact solution and the solution of MC\(_2(-0.21)\) compared to the solutions of MP and PS around the top of the soliton. (We omit the solution of DVD as it is the least accurate.) The approximate solutions have been reconstructed using cubic spline interpolation of the values at the grid points denoted with markers. This figure shows how well the solution of MC\(_2(-0.21)\) matches both the phase and the amplitude of the soliton.

Figure 3 shows the difference between the exact solution and the solutions of MC\(_4\)(0.06) and FD\(_4\) (with \(\varDelta t=0.5^2\)). The error in MC\(_4\) is roughly 30 times smaller, is located mainly at the peak of the soliton, and can be ascribed to a small phase error. The error of FD\(_4\) is more widespread.

Interaction of two solitons

We now study the interaction of two solitons over \(\varOmega =[-150,150]\times [0,50]\). The exact solution over \({\mathbb {R}}\) is [15],

where

Solution of EC\(_2(-0.18)\) on \(\varOmega \) and at the time \(T=50\)

Full size image

Comparison of different schemes around the top of the two solitons

Full size image

Interaction of two solitons: error of fourth-order methods

Full size image

We obtain the initial conditions from (46) setting

and solve this problem with \(\varDelta x=\varDelta t=0.5\). The optimal values of the free parameters for each of the families MC\(_2(\lambda _1)\), EC\(_2(\lambda _2)\), MC\(_4(\lambda _3)\), EC\(_4(\lambda _4)\) are \(\lambda _1=-0.19,\) \(\lambda _2=-0.18\), \(\lambda _3=0.06\), \(\lambda _4=0.02\). The results in Table 3 are consistent with those in Table 2, and analogous remarks apply.

Figure 4 shows the solution of the most accurate second-order scheme, EC\(_2(-0.18)\), on the whole domain \(\varOmega \) (upper plot) and at the final time (lower plot). The schemes MP and DVD produce oscillations (amplitude \(\simeq 0.005\)), where the exact solution is flat. These do not occur in the solution of EC\(_2(-0.18)\).

Figure 5 shows how the different schemes approximate the peak of the two solitons (omitting the least accurate solution of MP). The solution of EC\(_2(-0.18)\) best reproduces the speed and the amplitude of the two waves.

Finally, Fig. 6 compares the difference between the exact solution and the approximations given by MC\(_4\)(0.06) and FD\(_4\) (with \(\varDelta t=0.5^2\)). Just as for the single soliton, the error of FD\(_4\) has a higher amplitude and spreads far from the final location of the two solitons.

This section briefly demonstrates that the novel strategy described in Sect. 3.3 is practicable for PDEs with more than two independent variables. We seek to preserve two conservation laws,

of the pKP equation,

Note: Throughout this section, \(H_i\) and \({\widetilde{H}}_i\) are components of conservation laws, not Hamiltonians.

The characteristics \({\mathcal {Q}}_1=1, {\mathcal {Q}}_2=u_x,\) correspond, respectively, to the conservation laws with components

We introduce a uniform grid in space with nodes \((x_m,y_n)\) and use \(U_{m,n}(t)\) to denote a semidiscrete approximation of \(u(x_m,y_n,t)\). In this section, \(D_n\) and \(\mu _n\) are the forward difference and forward average operators acting on the second index, respectively.

The approach in Sect. 3.3 can be applied to a full 15-point rectangular stencil; this yields a wide range of families of methods. For brevity, we present here only those schemes for which all spatial derivatives are approximated on a one-dimensional spatial stencil consisting of three and five points, respectively, for the y- and x-derivatives. There are just two one-parameter families, both of the form

(so \(\widetilde{{\mathcal {Q}}}_1=1\)), that preserve semidiscrete versions of (48) and (49).

The first family is defined by

the semidiscrete version of (49) is

The second family has \({\widetilde{G}}_1\) and \({\widetilde{H}}_1\) as above, together with

The semidiscrete version of (49) is

For both families, \(\alpha =O(\varDelta x^2,\varDelta y^2).\) Let MC\(_1(\alpha )\) and MC\(_2(\alpha )\) denote the two families of fully discrete schemes obtained by applying implicit midpoint in time to (50) with (51) and (52), respectively.

Numerical test

As a brief test of the above schemes for the pKP equation, we use the following travelling wave solution of (47) [2]:

We apply methods MC\(_1(0)\) and MC\(_2(0)\) to the pKP equation on the domain \(\varOmega =[-0.5,0.5]\times [-10,10]\times [0,5]\) with initial and Dirichlet boundary conditions given by (53), using step lengths \(\varDelta x=0.01, \varDelta y=0.2\) and \(\varDelta t=0.05\).

Initial condition (left) and solution of MC\(_1(0)\) at time \(t=5\)

Full size image

Figure 7 shows the profile of the wave at the initial time and the numerical solution of MC\(_1(0)\) at the final time \(t=5\). Both MC\(_1(0)\) and MC\(_2(0)\) simulate the motion of the wave to the required accuracy, with a maximum absolute error in the solution of \(3.40\times 10^{-4}\) and \(3.41\times 10^{-4}\), respectively.

In this paper, we have introduced a new approach to constructing finite difference approximations to a system of PDEs that preserve multiple conservation laws. This is based on discretizing one dimension at a time, using semidiscrete Euler operators to find constraints that simplify the remaining symbolic computations. This is much cheaper than the approach introduced in [9] and can be iterated to apply to PDEs with more than two independent variables. We have proved that any symplectic Runge–Kutta method preserves local conservation laws with quadratic density and that the AVF method preserves the local conservation law of the energy under milder conditions than the skew-adjointness of the discrete operator \(\widetilde{{\mathcal {D}}}\).

The new strategy has been applied to obtain methods that preserve either the momentum or the energy of the Boussinesq equation. These are obtained as families that depend on a number of free parameters. Numerical tests have shown that the new schemes are competitive with respect to other methods in the literature and confirmed their conservation properties. Very accurate solutions can be obtained by selecting optimal parameter values. However, these values depend strongly on the choice of initial condition.

Finally, we have given an example that the new approach is practicable for PDEs with three independent variables, by finding two new families of schemes that preserve two conservation laws for the pKP equation.