A new technique for preserving conservation laws

This paper introduces a new symbolic-numeric strategy for finding semidiscretizations of a given PDE that preserve multiple local conservation laws. We prove that for one spatial dimension, various one-step time integrators from the literature preserve fully discrete local conservation laws whose densities are either quadratic or a Hamiltonian. The approach generalizes to time integrators with more steps and conservation laws of other kinds; higher-dimensional PDEs can be treated by iterating the new strategy. We use the Boussinesq equation as a benchmark and introduce new families of schemes of order two and four that preserve three conservation laws. We show that the new technique is practicable for PDEs with three dependent variables, introducing as an example new families of second-order schemes for the potential Kadomtsev-Petviashvili equation.


Introduction
Consider a system of q partial differential equations (PDEs), where A is a row vector, u has components u α , α = 1, . . . , q, and square brackets around a differentiable expression denote the expression and finitely many of its derivatives 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 Q such that Div F = AQ, (3) in which case Q is called the characteristic. The space of total divergences forms the kernel of the Euler operator, E, whose α-th entry is Hence if and only if there exists 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 [u] x denotes u and its spatial derivatives only, 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 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 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 higherdimensional spaces. We discuss this and give an illustration.
The paper is organised 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 Section 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 Section 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 Section 7.

Conservative space discretizations
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 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 centered difference operators are respectively, where I is the identity operator. The semidiscretizations of u α (x,t) are given by the column vector U ∈ R M q with the (m + α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 Div The following result is crucial for obtaining semidiscretizations that preserve conservation laws (see [18] and [14] for analogous results in the continuous and totally discrete setting, respectively).

Theorem 1
The kernel of the semidiscrete Euler operator E U , whose α-th entry is , is the space of semidiscrete divergences (10).
such that E U (L) = 0, and consider the derivative Integrating by parts yields ) whose precise expression is not of importance. Substituting this into (11) and integrating over ε ∈ [0, 1] shows that L is a semidiscrete divergence. If L is of the form (10), E U (L) = 0 follows from the linearity of the Euler operator and from the fact that for any k (see, e.g., [12]), 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 Q ℓ , as follows: 1. Select a stencil that is large enough to support generic semidiscretizations for A and every Q ℓ , having the desired order of accuracy. These approximations depend on a number of free parameters to be determined.
2. Find some of the parameters by imposing consistency, up to the desired order of accuracy, p.
3. Use symbolic algebra to determine the values of the free parameters that satisfy for ℓ = 1. This guarantees that the first conservation law is locally preserved. As both A and Q ℓ are accurate to order p, the discrete conservation law has the same order of accuracy.
4. Iterate the previous step, replacing Q 1 with Q ℓ , to obtain further constraints on the parameters. If (12) has no solution for some ℓ, 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.

Time integration
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 and let u m,n ∈ R q be the column vector with entries u α m,n ≈ u α (x m , t n ), α = 1, . . . , q.

Conservation laws with quadratic density
Here attention is restricted to PDEs of the form where g is linear homogeneous in [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 (x) be an invertible operator such that and h(x, t, U) are two column vectors whose (m + αM − B)-th entry is a spatial discretization of the α-th component of g( 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 Section 2. The flux and density of (16) have the general form 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 u n = P (x) g n . Moreover, Using (19) to eliminate g n from the first sum and rearranging, gives The condition of symplecticity, and (17) give which is an approximation of (14).
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 [u] x .

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): 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 ([u] x ) on such a stencil by H(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([u] x ) as H(µ m 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 D is a skew-adjoint approximation of 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 respectively.

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.

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. Let P = r i=1 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 U or U t . Assuming that the stencil has N points in the time dimension, discretize P using where i j is the i-th 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 then 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 α j only satisfy the requirements for consistency. For example, in a one-step method linear quantities are approximated by and quadratic quantities by 2. The values of the parameters α j are obtained by solving where A and Q ℓ are the approximations of the PDE and the characteristic obtained after steps 1 and 2, and E u is the difference Euler operator, 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: • M + r − 1 r 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 semidscretizing first, the strategy in [9] introduces NM + r − 1 r 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 k-th 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.

The Boussinesq equation
We consider here the (Good) Boussinesq equation cast as a system of two PDEs This system can be written in Hamiltonian form, System (31) has infinitely many independent conservation laws [1]. The first four are D 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.

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 Q 1 = Q 1 and Q 2 = Q 2 . The linear and quadratic terms in (31) and (36) are approximated as and the coefficients α i and β 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 Section 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 Q 3/4 and density G 3/4 of the remaining preserved conservation law. The corresponding flux 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 (α, β, γ, ξ) are O(∆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 ∈ {0, 1/3, 1/2}.
The first three families are Three remaining are given by F 1 , F 2 and G 2 as in (39) together with . In the numerical tests section we limit our investigations to the semidiscretions obtained from (38) with (39) and s = α = β = ξ = 0, γ = λ 1 ∆x 2 ; we use MC 2 (λ 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 Section 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, F 1 and G 2 are defined as in (39) and

The resulting system of ODEs can be written in the form
with The operator 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 (λ 2 ) to denote the schemes with α = γ = 0 and β = λ 2 ∆x 2 . Again, the iterative technique from Section 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 (α, β, γ, ξ) that are all O(∆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 F j are second-order accurate, but D m 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 F 1 , F 2 and G 2 as in (42), together with We use MC 4 (λ 3 ) to denote the schemes obtained by applying the Gauss-Legendre method of order four to (42), with α = β = γ = 0 and ξ = λ 3 ∆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 ∈ N, let and define the operators ν ± m = I ± ∆x 2 6 D 2 m S −1 m and the functions The family of semidiscretizations and conservation laws is as follows: The systems of ODEs defined by (38) with (43) amounts to with The operator D, although not skew-adjoint, satisfies (23). We use EC 4 (λ 4 ) to denote the family of schemes obtained by applying the AVF method of order four (see [19]) to (44) with α = λ 4 ∆x 4 and β = γ = 0.

Numerical Tests
In this section we solve a couple of benchmark problems to show the effectiveness and conservation properties of the numerical methods in Section 4.
The results are compared with the following second-order structure-preserving methods: • The multisymplectic scheme for (30), [22] and equivalent to the well-known Preissmann scheme.
• The energy-conserving scheme for (31) in [16], , 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 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 Section 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 ∆t 2 . We consider (x, t) ∈ Ω ≡ [a, b] × [0, T ] and periodic boundary conditions. We introduce on Ω 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 ≃ u(a + i∆x, j∆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 · denotes the Euclidean norm. We also compare the errors in the global invariants (37) defined as where α = 1, 2, 3, 4. For the methods introduced in this paper, G α is given in Section 4. For all the other schemes, we set

Single soliton
For the first problem we set Ω = [−60, 60] × [0, 25] and the initial conditions given by the single soliton solution over R, We first examine the convergence and stability of the schemes found in Section 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 ∆x = ∆t = k/10, for k = 1, 2, . . . , 9.
The results in Table 1 and Figure 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 stepsizes. Table 2 shows the error in the conservation laws and the solution for the different methods with ∆x = ∆t = 0.5. For this problem, the values of the free parameters that minimize the 6.21e-15 1.57e-12 1.47e-08 5.57e-14 1.94e-04 FD 4 (∆t = 0.5 2 ) 9.97e-12 0.0036 1.07e-04 2.84e-04 0.0038 error in the solution of MC 2 (λ 1 ), EC 2 (λ 2 ), MC 4 (λ 3 ), EC 4 (λ 4 ) are λ 1 = −0.21, λ 2 = −0.20, λ 3 = 0.06, λ 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 optimising 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 timestep ∆t = 0.5 2 , even for non-optimal parameters.
In Figure 2 Figure 3 shows the difference between the exact solution and the solutions of MC 4 (0.06) and FD 4 (with ∆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 Ω = [−150, 150] ×[0, 50]. The exact solution over R is [15], where We obtain the initial conditions from (46) setting and solve this problem with ∆x = ∆t = 0.5. The optimal values of the free parameters for each of the families MC 2 (λ 1 ), EC 2 (λ 2 ), MC 4 (λ 3 ), EC 4 (λ 4 ) are λ 1 = −0.19, λ 2 = −0.18, λ 3 = 0.06, λ 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 Ω (upper plot) and at the final time (lower plot). The schemes MP and DVD produce oscillations (amplitude ≃ 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, Figure 6 compares the difference between the exact solution and the approximations given by MC 4 (0.06) and FD 4 (with ∆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.

The potential Kadomtsev--Petviashvili (pKP) equation
This section briefly demonstrates that the novel strategy described in Section 3.3 is practicable for PDEs with more than two independent variables. We seek to preserve two conservation laws, of the pKP equation, u xt + 3 2 u x u xx + 1 4 u xxxx + 3 4 u yy = 0.
Note: Throughout this section, H i and H i are components of conservation laws, not Hamiltonians.
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 µ n are the forward difference and forward average operators acting on the second index, respectively.
The approach in Section 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 D m F 1 + D n G 1 + D t H 1 = 0 (50) (so Q 1 = 1), that preserve semidiscrete versions of (48) and (49). The first family is defined by Figure 7: Initial condition (left) and solution of MC 1 (0) at time t = 5.

Numerical test
As a brief test of the above schemes for the pKP equation, we use the following travelling wave solution of (47) [2]: u(x, y, t) = 2 tanh(x + y − 7 4 t + 5) + 2.  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 × 10 −4 and 3.41 × 10 −4 , respectively.

Conclusions
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 discretising 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 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.