{'fraction_non_alphanumeric': 0.04564032697547684, 'fraction_numerical': 0.021984146643547187, 'mean_word_length': 4.27039817232376, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On structure-preserving discontinuous Galerkin methods for Hamiltonian partial differential equations: Energy conservation and multi-symplecticity

In this paper, we present and study discontinuous Galerkin (DG) methods for one-dimensional multi-symplectic Hamiltonian partial differential equations. We particularly focus on semi-discrete schemes with spatial discretization only, and show that the proposed DG methods can simultaneously preserve the multi-symplectic structure and energy conservation with a general class of numerical fluxes, which includes the well-known central and alternating fluxes. Applications to the wave equation, the Benjamin-Bona-Mahony equation, the Camassa-Holm equation, the Korteweg-de Vries equation and the nonlinear Schr\"odinger equation are discussed. Some numerical results are provided to demonstrate the accuracy and long time behavior of the proposed methods. Numerically, we observe that certain choices of numerical fluxes in the discussed class may help achieve better accuracy compared with the commonly used ones including the central fluxes.


Introduction
Structure-preserving numerical methods, which preserve certain structures and invariants of the model in the discrete level, have gained many attention in the simulation of mathematical models. It is well-known that such schemes tend to achieve better long time behavior in terms of stability and accuracy for Hamiltonian dynamics [23]. For the multisymplectic HPDEs (1.1), numerical methods preserving (1.3) are referred as multi-symplectic integrators. There have been many studies on such methods, with examples including the Preissman box schemes [3,46,1,17], Euler box schemes [32,17], diamond schemes [31], spectral methods [4,8], multi-symplectic (partitioned) Runge-Kutta (RK) methods [35,21,36] and recently, DG methods [40,6]. These methods have been successfully applied to various equations including the Hamiltonian wave equation [3,31], the BBM equation [37,24], the CH equation [17], the KdV equation [46,1], the Schrödinger equation [8,40] and the Dirac equation [20]. Recently, there have been increasing interests in designing local energy conserving numerical schemes for the continuous dynamical systems. In the case that the associated energy and momentum functionals of (1.1) are quadratic, some of these multisymplectic methods also preserve the local conservation of energy (1.4) and momentum (1.5) at the discrete level, while in general this is not the case. Usually, particular discretizations have to be constructed to preserve the energy or the momentum of (1.1), such as averaged vector fields (AVF) methods [30,33,7] and invariant energy quantization (IEQ) methods [5].
In this paper, we investigate the structure-preserving property of the DG spatial discretization for the multi-symplectic HPDEs (1.1). The DG method is a class of finite element methods using discontinuous piecewise polynomial spaces. It was first introduced by Reed and Hill [34] for solving the transport equation and then received its major development in a series of work by Cockburn et al. in [14,13,12,11,16] for hyperbolic conservation laws. For equations containing high order spatial derivatives, the local discontinuous Galerkin (LDG) methods were proposed by Cockburn and Shu in [15]. The DG method also finds its strength on preserving structures of the continuum equations (sometimes with suitable limiters), such as the positivity or other physical bounds [45], the hydrostatic balance [42,43], the entropy inequality [9,38] and the asymptotic limits [18]. Recently, there have been many studies in designing DG and LDG methods which can numerically preserve the energy or Hamiltonian structure of the model. Energy conserving LDG methods have been designed for the generalized KdV equation [2,29], the acoustic wave equation [10], the CH equation [28], the Degasperis-Procesi equation [22], the nonlinear Schrödinger equation [27], the improved Boussinesq equation [25] and so on.
This paper can be considered as a step further along this track, and our goal is to design structure-preserving DG methods for the multi-symplectic HPDEs (1.1) which can preserve both the multi-symplectic structure and associated local energy conservation. This work is also a reinterpretation and generalization of the earlier framework by Tang et al. [40] on DG methods with alternating fluxes for HPDEs. When K has a blocked structure, it was shown in [40] that DG methods for (1.1) with alternating fluxes can be rephrased as partitioned RK methods. Then with suitable quadrature rules and time integrators, the fully discrete methods can be interpreted as space-time partitioned RK methods, whose multisymplecticity follows from [36]. In this paper, our attention is particularly on method-of-lines DG schemes. We derive the multi-symplecticity directly from the weak formulation. The analysis covers a very general class of numerical fluxes including, but not limited to, the well-known central and alternating fluxes studied in [6]. We observe that, in the case when neither central nor alternating fluxes give optimal convergence rate, the proposed wider range of numerical fluxes selections may recover the optimal convergence numerically. Furthermore, these DG methods are shown in [6] to be energy conserving for two particular examples of HPDEs, namely the Schrödinger equation and KdV equation. Motivated by this, we are able to prove that semi-discrete DG methods simultaneously conserve the associated energy for general multi-symplectic HPDEs (1.1). This result also indicates a rather general approach of designing DG schemes that conserves a certain Hamiltonian invariant of HPDEs. If one can construct a multi-syplectic system with the aimed invariant as the associated energy, it will be automatically conserved by the proposed DG approximation applied to such multisymplectic formulation. For example, we can show that the proposed multi-symplectic DG scheme, with central fluxes, for the KdV equation retrieves the Hamiltonian preserving DG scheme presented in [29].
For applications, we especially consider the Hamiltonian wave equation, the BBM equation and the CH equation in this paper. The choices of numerical fluxes and the corresponding implementation procedures are discussed in details. Since the focus of the paper is on spatial discretization, for numerical tests, we usually apply high order RK time integrators with small time steps to reduce temporal error, so that fully discrete schemes faithfully approximates method-of-lines schemes. In principle, one can apply a suitable symplecitic RK method to preserve multi-symplecticity, or an AVF method to preserve the energy. However, preserving both properties simultaneously can be difficult for problems with non-quadratic energy functionals. In the numerical tests, we also observe certain choices of numerical fluxes may help improve the accuracy while still conserving its corresponding energy. Other tests, such as multi-wave interactions, are also provided to illustrate the performance of the DG schemes.
The rest of the paper is organized as follows, in Section 2 we state DG methods for multi-symplectic system (1.1). The preservation of multi-symplectic structure and the local energy conservation are proved in Section 3. After that, application to various HPDEs and the implementation procedure are discussed in Section 4, and numerical tests are provided in Section 5. Finally, we close the paper with conclusions in Section 6.

The DG scheme
Consider a quasi-uniform partition of the spatial domain Ω ), and the mesh size is denoted by with h = max h j for j = 1, 2, ..., N being the maximal mesh size. We use P k (I j ) to represent the linear space spanned by polynomials of degrees no more than k on I j . Let to represent the function limit from the right and left respectively. We denote the average and jump of the functions at the cell interfaces by . The semi-discrete DG approximation to (1.1) is given as follows: Seek the numerical solutions z z z h ∈ V V V h , such that holds for all test functions ϕ ϕ ϕ ∈ V V V h . The hatted terms, Kz z z h , are the numerical fluxes defined on the element interfaces, and are the key component in designing the DG methods. In this paper, we choose the family of the numerical flux Kz z z h to be for any m × m real symmetric matrix A and real anti-symmetric matrix B. Some comments on the choice of the numerical flux Kz z z h are given below. We will show that this family of numerical flux includes the central flux and the alternating flux, wellknown in the LDG methods applied to equations with high order derivatives. Since K is anti-symmetric, there exists an orthogonal matrix Q such that with α ∈ [−1/2, 1/2], the numerical flux in (2.2) reduces to which retrieves alternating fluxes with α = ± 1 2 , and central fluxes with α = 0. Second, from the practical point of view, B = 0 should work for most situations. While in numerical tests, we notice that a nonzero B may help improve the accuracy when A = 0. We therefore include B in our analysis for completion. Further comments on choices of numerical fluxes are postponed to Remark 6.1 in the conclusion.
The DG scheme (2.1) for HPDEs (1.1) also relates to an LDG method for the associated scalar equation. The DG discretization is applied to a first order system reformulated from the original scalar equation and all auxiliary variables will be eliminated. In our cases, the first order system is particularly given as the multi-symplectic system and the elimination procedure will be detailed in Section 4.

Properties of the DG scheme
In this section, some properties of the proposed semi-discrete DG scheme (2.1) are investigated. More specifically, we are interested in demonstrating that the multi-sympleciticy and the energy conservation are both preserved by our methods. We focus on the spatial discretization only, and will comment on the appropriate temporal discretization at the end of this section.

Multi-symplecticity
Applying the exterior derivative to the DG scheme (2.1) yields the following variational equation for the one-forms of z z z h (3.1) To facilitate our discussion, the following equalities are provided.
Proof. We provide the proof of (3.2) only, and skip that of (3.3) which follows along the similar line. Introduce the notation of and it yields that by combining corresponding terms in the subtraction. Using the definition of Kz z z h in (2.2), one can obtain It can be verified that since the matrices K, B are anti-symmetric and A is symmetric. After substituting (3.6) into (3.5), one can then obtain (3.2).
Suppose dz z z h , dz z z h ∈ V V V h both satisfy (3.1), then we have the following theorem.
Theorem 3.1 (Conservation of multi-symplecticity). For any one-forms dz z z h , dz z z h ∈ V V V h satisfying the variational equation (3.1), with fluxes defined in (2.2), we have the semidiscrete version of the multi-symplectic conservation laws where F is defined in (3.4) and Proof. Using the anti-symmetric property of the matrix M and the variational equation (3.1), it can be shown that With the anti-symmetry of K and the symmetry of ∇ z z zz z z S(z z z h ), we have With fluxes defined in (2.2), one can then apply Lemma 3.1 to complete the proof.
Remark 3.1. When the anti-symmetric matrix B is chosen as zero, ω h,j defined in (3.7) reduces to I j M dz z z h · dz z z h dx. Also, when z z z h ,z z z h are continuous across cell interfaces, we have Remark 3.2. As can be seen from the proof, other choices of fluxes besides (2.2) also preserve the mulit-symplecticity, as long as is equal to the time derivative of a term consistent with zero.

Local energy conservation
In this subsection, the local energy conservation property of the proposed DG scheme (2.1) is explored. Recall that the continuous energy E(z z z) = S(z z z) − 1 2 Kz z z x · z z z is defined in (1.4), and we have the following theorem on the local conservation of the discrete energy.
Theorem 3.2 (Conservation of energy). The numerical solution z z z h of the semi-discrete DG scheme (2.1) satisfies the local energy conservation law in the form of As a consequence, it conserves the total energy Proof. Following the definition of the energy E(z z z), we have (3.10) The numerical solution z z z h satisfies the variational equation (2.1). Taking the test function ϕ ϕ ϕ = (z z z h ) t leads to where the last equality follows from the anti-symmetry of M and K. By substituting (3.11) into (3.10), we have On the other hand, using the product rule, it can be shown that (3.12) where the last equality follows from Lemma 3.1. The local energy conservation property (3.8) can be obtained after rearranging terms of (3.12). The global energy conservation property (3.9) follows from summing it over all elements and utilizing the definition of the numerical flux in (2.2) and the periodic boundary condition.
The global energy term E h in (3.9) can be further simplified for some special choices of the matrices A and B in the definition of the numerical flux, which is summarized in the following corollary.
if m is even. (3.14) In particular, for A defined in Proof. Without loss of generality, we assume m is odd in this proof. We first derive the variational form associated with the last m 2 equations. Note the DG scheme (2.1) has the following global form after summing over all the elements: Recall that B = βM , which leads to B[z z z h ] t · [ϕ ϕ ϕ] = 0 following the steps in (3.16). Similar computation as in (3.17) yields Moreover, since We now turn to simplify the energy functional E h defined in (3.9). One can apply integration by parts to get On the other hand, it can be shown that Therefore, we can rewrite the energy functional E h as (3.20) One can easily verify that (3.22) The equality (3.13) is obtained by substituting (3.21) and (3.22) into (3.20), and using the definition of A in (3.14). This completes the proof.

Application and implementation
In this section, we present the DG discretization of a few particular examples of the multisymplectic HPDEs, including the Hamiltonian wave equation, the BBM equation, the KdV equation, the Schrödinger equation and the CH equation. Our main attention is on various numerical fluxes and discrete energy conservation indicated by Theorem 3.2. We detail the choices of numerical fluxes and implementation of our DG method for the Hamiltonian wave equation, the BBM equation and the CH equation. It will be shown that for these models, the proposed methods with appropriate numerical fluxes recover some existing energy conserving DG methods studied in the literature. We have also listed the application to the Schrödinger equation and the KdV equation, while these schemes (with central and alternating fluxes) have been studied in [6].
For all of these examples of HPDEs, the original equation involves only one unknown u, and they can be rewritten in the form of HPDEs by introducing auxiliary variables z z z. The DG method (2.1) is presented with the unknown variable z z z h . In the numerical implementation of the DG scheme, we want to eliminate all the intermediate unknowns and convert the system back into an equation involving only one unknown u h , which will be explained in this section. The following operators are useful in this procedure to simplify the notations.
• We also use Π : to represent the standard L 2 projection to the piecewise polynomial space V h .
With a given set of basis functions, matrices for the local operator D α and L can be assembled explicitly in each interval I j . The projection Π can also be implemented as a computer subroutine. To implement a DG scheme, for example (4.9) for the Hamiltonian wave equation, one simply replaces trial functions and operators with corresponding vectors and matrices or computer subroutines.

Hamiltonian wave equation
The Hamiltonian wave equation takes the form of and includes the linear second order wave, Sine-Gordon and Klein-Gordon equations as special cases. By introducing two auxiliary variables v and w, it can be formulated as a multi-symplectic system (1.1) with z z z = (u v w) T and which can be expressed as By utilizing the fact that u x = w, the energy functional E of the continuum equation can be simplified as Applying DG discretization (3.1) to (4.5) gives the following scheme: The energy conserving LDG method (with alternating fluxes) for the second order wave equation has been extensively studied in [41], where the optimal error estimate, energy conserving, superconvergence properties are carefully analyzed. It was also demonstrated numerically, that the shape of the solution, after long time integration, is well preserved due to the energy conserving property. The DG method (4.6), with alternating fluxes α 11 = α 33 = 0, α 13 = ± 1 2 and β = 0 in (4.7), retrieves the same energy conserving LDG method in [41], which means that the proposed semi-discrete method is both multi-symplectic and energy conserving.
which can be further simplified as

Nonlinear KdV equation
Next, we consider the nonlinear KdV equation where η, ε are given parameters. It can be written as a multi-symplectic system (1.1) with which can be expressed as (4.11) The continuum equation is associated with the energy functional To be more specific, we look for u h , v h , w h ∈ V h , such that for all ϕ i ∈ V h , i = 1, 2, 3, 4, we If we choose α 1 = 0, and α 2 = 0 (leading to central flux) or α 2 = ± 1 2 (leading to alternating flux), the proposed methods reduce to the same ones that was studied in [40].
we have the following discrete energy conservation from Corollary 3.1.

Proposition 4.2.
For any α 1 and α 2 , the DG scheme (4.13) with numerical fluxes (4.14) conserves the discrete energy As for the implementation, we look into the DG scheme in the operator form.
When α 1 = 0, one can eliminate all auxiliary variables to obtain The "energy conserving" DG method, which conserves the Ω u 2 h dx (the momentum in the multi-symplectic community), have been recently proposed and studied in [2]. Suboptimal error estimate was provided. Later, the Hamiltonian preserving DG method, which conserves the energy (4.12), was presented in [29], and it was shown in the paper that the Hamiltonian preserving DG methods have a slightly improved long time behavior over those in [2]. The DG method (4.13), with central fluxes, retrieves the same Hamiltonian conserving DG method in [29], which means that semi-discrete method is both multi-symplectic and energy conserving. To see connections between the scheme (4.13) and that in [29], we apply 1 2 ∂ t to (4.13d), take ϕ 4 = ϕ 1 and then combine it with (4.13a) to get (4.15) and (4.13b). Then (4.15), (4.13b) and (4.13c) recover the scheme in [29].

BBM equation
The BBM equation with σ being a given positive parameter, models long shallow water waves of small amplitude and is widely studied. It can be reformulated as (by introducing V (u) = −u 3 /6) Then the BBM equation (4.16) can be expressed as After simplification, the corresponding energy functional is Applying DG discretization (3.1) to (4.16) yields the following scheme: Find φ h , u h , v h , w h , p h ∈ V h , such that for all ϕ i ∈ V h , i = 1, 2, 3, 4, 5, we have take the form of . One can apply Corollary 3.1 to show the proposed DG method conserves a modified energy, given in Proposition 4.3.
For implementation, the scheme (4.17) and (4.18) can be rewritten as the following system One can eliminate v h , w h and p h to obtain Note although α 2 appears in the weak formulation, it eventually gets cancelled and does not play a role in the numerical scheme. Case 1 (Central flux): Suppose α 0 = α 1 = 0, we can eliminate φ h to get Numerically, we see that this scheme is (k + 1)th order accurate for even k on uniform mesh, but is only k-th order for odd k or on nonuniform mesh. Case 2 (Generalized alternating flux): Suppose α 0 = 0 and α 1 = 0. We solve for φ h from u h in (4.19e) for further reducing the scheme. The method can then be formally written as However, the kernel of D α 1 is one dimension and D α 1 −1 is only well-defined on {v h ∈ V h : Ω v h dx = 0}. Therefore, when implementing the numerical scheme, we do time marching for u h − Ω u h dx. We drop the last cell average and perform the matrix inversion on a smaller vector in R N (k+1)−1 , and finally retrieve the last cell average using the zero average condition. To have the other matrix inverse well-defined, we need the polynomial degree k to be even and the number of spatial mesh cells N to be odd. Although the implementation of this case can be quite involved, numerically it still suffers order reduction as in Case 1.
Case 3: We then consider α 0 = 0 and α 1 = 0. Then we can obtain This scheme seems to retrieve optimal convergence rate for u h numerically.

CH equation
The nonlinear CH equation, a bi-Hamiltonian model for waves in the shallow water, takes the form of u t − u xxt + 3uu x − 2u x u xx − uu xxx = 0, and corresponds to the multi-symplectic system (1.1) with z z z = (u φ ρ v w) T , and The CH equation can then be expressed as The associated continuous energy functional is given by The variational form of the DG scheme (3.1) for the multi-symplectic CH equation is given as follows: Find u h , φ h , w h , v h , ρ h ∈ V h , such that for all test functions ϕ i ∈ V h , i = 1, 2, 3, 4, 5, we have One can apply Corollary 3.1 to show the proposed DG method conserves a modified energy, as explained below.
For implementation, the DG scheme can be written in the following operator form.

Nonlinear Schrödinger equation
Consider the nonlinear Schrödinger equation where i is the imaginary unit and α is a given positive parameter. Let u = p + iq. Then

This continuum equation preserves the energy functional
We omit the detailed schemes in this section to save space. With alternating fluxes, the DG method (3.1) applied to (4.22), is the same as the original LDG method for the NLS equation introduced in [44]. Later, the same DG method, with central or alternating fluxes (3.1), has been studied in [6] and [40] for its multi-symplectic property. The energy conservation has also been investigated in these papers, which can also be obtained through Corollary 3.1 as below.
As has been pointed out in [6], this DG scheme also conserves the total charge C h = Ω p 2 h + q 2 h dx. This conservation property relies on special structures of (4.22), and is hence not covered under the current framework.

BBM-KdV equation
After simplification, the associated energy functional is given by To avoid complications in implementation, we consider the DG method with central fluxes for (4.23). The scheme actually retrieves that for KdV equation (4.10) with σ = 0, V (u) = − η 6 u 3 and that for BBM equation (4.16) with ν = 0. Furthermore, we have the following energy conservation property.

Numerical tests
In this section, we provide some numerical results to demonstrate the behavior of the proposed multi-symplectic and energy conserving DG methods for the Hamiltonian wave equation, the BBM equation and the CH equation. We refer to [6] for the performance of the DG methods for the KdV and Schrödinger equations, when central or alternating fluxes are used. In the accuracy tests, both uniform and nonuniform meshes are considered. For nonuniform meshes in all tests, the mesh size ratios are set as 2 : 1 : 2 : 1 : · · · , i.e., we have ∆x 2j−1 = 2∆x 2j , j = 1, 2, · · · N 2 with N being the number of grid points. For simplicity, periodic conditions are used for all tests. Various RK temporal discretizations are used in the tests, with details provided in each test.

Hamiltonian wave equation
The DG method with alternating numerical fluxes for the second order wave equation has been studied in [41], and here we provide some additional numerical results when other choices of numerical fluxes are considered.
Example 5.1.1 (Accuracy test). In this example, we examine the accuracy of the DG scheme (4.6) for the wave equation (4.4) with different choices of the numerical fluxes (4.7). Consider the setup of the problem as u(x, 0) = e sin(x) , u t (x, 0) = cos(x)e sin(x) .
The exact solution to the problem is u(x, t) = e sin(x+t) . For DG scheme with P k elements, we apply (k + 1)th order RK method to (4.8) for time discretization. The time step is taken as ∆t = 0.01∆x. To reveal the convergence rate, we carefully choose the initial condition u h , such that w h also has the optimal convergence rate. We start with examining the accuracy with central fluxes (α 11 , α 13 , α 33 , β) = (0, 0, 0, 0). From Table 5.1, one can see that, for uniform meshes and for both u h and w h , the scheme is kth order accurate if k is odd and (k + 1)th order if k is even. However, for nonuniform meshes, the scheme degenerates to kth order for all cases. These observations are consistent with our understanding when central fluxes are used.
Then we study the effect of numerical fluxes on accuracy. Only P 1 and P 2 schemes are presented for simplicity. Tests in Table 5.2 indicate that the accuracy of u h can be improved if any of the following hold: α 11 ≥ 0, α 13 = 0 or β = 0; while optimal rates of w h are retrieved when α 33 < 0 or α 13 = 0.
Example 5.1.2 (Error and energy change). In this test, we focus on the numerical error and energy change (∆E h (t) := E h (t)−E h (0)) after long time simulation. We consider the problem (5.1a) together with the initial condition u(x, 0) = sin(cos(x)) and u t (x, 0) = 0. The exact solution of the problem is u(x, t) = 1 2 (sin(cos(x + t)) + sin(cos(x − t))). P 3 elements with 100 mesh cells are used for numerical simulations. We use the fifth order RK method for time marching and the time step is set as ∆t = 0.01∆x to reduce the temporal error. The final time is set as T = 200π, which corresponds to 100 period. The plots for L 2 error and Uniform mesh  the energy change are given in Figure 5.1. With central fluxes, the scheme has numerical error around 10 −7 , and the error is around 10 −8 for other numerical fluxes. It can also be seen that the energy change remains at a small magnitude of around 10 −13 ∼ 10 −12 for all these tests, and the minor energy change seems to be due to time discretization which is not energy conserving.
We start with the accuracy test. Again, for P k elements, explicit RK method of order k + 1 is used for time integration. The time step is taken as ∆t = 0.5∆x and we compute to T = 1. From Table 5.3, we observe a similar convergence pattern as that in Table 5.1 if α 0 = 0 (which includes central fluxes). While if we set α 0 = 1 4 and α 1 = α 2 = 0 for numerical fluxes, it achieves the (k + 1)th order convergence rate for u. Other choices of α 0 > 0 seem to retrieve the optimal rate as well. We remark that a large α 0 may result in a more restricted time step, and a small α 0 may require a more refined mesh to observe the actual convergence rate. When α 0 < 0 is used, the method seems to be unstable numerically.

Uniform mesh
Nonuniform mesh  Example 5.2.2 (Energy conserving property). We use the same test problem as in the previous example, to compare the performance of E 3 -conserving scheme (which preserves E 3 = 1 6 Ω u 3 dx) in this paper and the E 2 -conserving scheme [26] (which preserves E 2 = Ω (u 2 + σu 2 x )dx) for long time simulation. Here E 3 and E 2 denote the third and second invariant of the BBM equation, respectively. P 2 polynomials with 10 mesh cells (∆x ≈ 9.2237 × 10 −2 ) are used in the test. To reduce the temporal error, we apply the fifth order RK method with a small time step ∆t = 0.01∆x for time discretization. The L 2 error, ∆E h,3 and ∆E h,2 are plotted in Figure 5.2, and the snapshots at T = 200, 1000, 3000, 5000 are given in Figure 5.3. By comparing the Figure  5.2(a) and Figure 5.2(d), one can tell the L 2 error of the E 3 -conserving scheme grows slower than that of the E 2 -conserving scheme. According to Figure 5.3, although both methods preserve the profile of the traveling wave well after very long time, the E 3 -conserving scheme commits slightly smaller phase error compared with the E 2 -conserving scheme.
In this example, we perform numerical simulations for the soliton and multi-soliton interaction of the BBM equation. The initial data are given in Table 5.4. We fix σ = 11 100 2 and use P 4 elements with central fluxes for spatial discretization. For time integration, the fifth order RK method with ∆t = 0.05∆x is used. Settings of the numerical tests and the corresponding results are documented in Table 5

CH equation
Example 5.3.1 (Accuracy test). We start with the accuracy test with the fabricated solution. The exact solution is set as u(x, t) = sin(x + t) and the source term is computed accordingly. We compute to T = 1 with the time step ∆t = 0.01∆x using (k + 1)th order RK scheme. Numerical fluxes with α 0 = 0 and α 0 = 3 are tested under this setting. The numerical error is documented in Table 5.5, where the same even-odd phenomenon is observed when α 0 = 0 and the optimal convergence rate for u h seems to be retrieved when α 0 is chosen to be 3.

Uniform mesh Nonuniform mesh
In the following numerical tests, we simulate the single traveling peakon, two-peakon interaction, three-peakon interaction and peakon-antipeakon interaction using multi-symplectic DG scheme with central fluxes. Except for the single traveling peakon, other test problems are taken from [28]. P 4 elements with ∆x = 0.075 are used to resolve the solution. The spatial domain is set as Ω = (0, 30). To ensure stability, we use the third order strongstability-preserving RK method with superviscosity stabilization after each time stage. See [39] for further studies on this stabilization approach. The time step is set as ∆t = 0.01∆x. Other settings of the tests are documented in Table 5.6. The profiles of numerical solutions are consistent with those in [28].

Conclusions
In this paper, we study the semi-discrete DG discretization of multi-symplectic systems, along with its application to various multi-symplectic HPDEs. With a general class of numerical fluxes, the semi-discrete DG schemes are shown to preserve the multi-symplectic structure and the local energy conservation simultaneously. For fully discrete schemes, one of (a) Solution profile.       the properties can also be achieved with suitable time integrators chosen. For applications, we particularly consider the wave equation, the BBM equation and the CH equation. The preserved energy functionals and their implementation methods are discussed. Numerically, we observe that different choices of numerical fluxes do have an influence on the accuracy of the schemes, and optimal convergence rate could be achieved with carefully chosen numerical fluxes. The DG scheme for each equation preserves its corresponding energy functional well in long time simulation with accurate time discretization, and enjoy the benefit of energy conserving methods like an improved long time behavior. Numerical simulations of multiwave interactions are also provided to illustrate the performance of the scheme. A general class of numerical fluxes has been discussed in the paper, and we would like to provide the following remark as a guideline on the choice of numerical fluxes.
Remark 6.1 (Choices of numerical fluxes). Although in principle any numerical fluxes of the form (2.2) preserve the multi-symplectic structure and the energy conservation, in practice, we recommend fluxes that are easy to implement and retrieve optimal convergence rates. Central fluxes, with A = B = 0, usually lead to a scheme that is easy to implement. While it may suffer order degeneration for odd order polynomials or on nonuniform meshes. We recommend to use alternating fluxes if the implementation is not an issue, since the convergence rates are usually optimal for both u h and its discrete derivative. While alternating fluxes may not work as desired for the BBM and CH equations. In this situation, we suggest to tune A such that the numerical flux for u h is no longer central, and then supplement extra fluxes to ensure symmetry of A. With this approach, the accuracy can be improved and the implementation is usually not that complicated.