High-Order Energy and Linear Momentum Conserving Methods for the Klein-Gordon Equation

: The Klein-Gordon equation is a model for free particle wave function in relativistic quantum mechanics. Many numerical methods have been proposed to solve the Klein-Gordon equation. However, efﬁcient high-order numerical methods that preserve energy and linear momentum of the equation have not been considered. In this paper, we propose high-order numerical methods to solve the Klein-Gordon equation, present the energy and linear momentum conservation properties of our numerical schemes, and show the optimal error estimates and superconvergence property. We also verify the performance of our numerical schemes by some numerical examples.


Introduction
The Klein-Gordon equation is a widely used model in relativistic quantum mechanics.It is used to describe the free particle wave function.Many numerical methods, including finite difference methods [1][2][3][4], finite element methods [5] and spectral methods [6,7], have been proposed to solve the Klein-Gordon equation.In particular, the authors in [2] proposed a simple second order centered difference in time and space, an explicit scheme to solve the Klein-Gordon equation.The finite difference schemes proposed in [1,4] used the same treatment for second order spatial and temporal derivatives, but different treatment for the source term of the equation.In [3], a finite difference method along with an operator splitting method was proposed to solve the equation on an unbounded domain.All of these finite difference schemes have second-order convergence in space and time.Some types of spectral methods were proposed in order obtain a higher order of convergence in space.In [6], a pseudo-spectral method was proposed for the spatial discretization, and Crank-Nicolson or leap-frog method was used for the temporal discretization.In [7], Legendre wavelets incorporated with a spectral method were used to solve Klein-Gordon and Sine-Gordon equations.As far as the family of finite element methods is concerned, Ref. [5] proposed Galerkin finite element methods for the Klein-Gordon equation with homogeneous Dirichlet boundary conditions.However, none of the aforementioned work focused on the conservation properties of energy and momentum.The advantage of high-order energy and linear momentum conserving methods is that they can lead to accurate and physically relevant solutions with relatively coarse mesh.Moreover, the structure-preserving property is a significant factor when judging the performance of numerical schemes [8].
In this paper, we propose high-order local discontinuous Galerkin (LDG) methods that have provable properties of energy and linear momentum conservation, optimal convergence and superconvergence.Discontinuous Galerkin (DG) methods belong to the family of finite element methods, and have drawn great attention since the work by Cockburn and Shu in [9][10][11][12].LDG methods, as a special type of DG method, were then proposed to solve diffusion equations [13], partial differential equations with high spatial derivatives [14,15], and wave equations [16][17][18][19].In this work, we propose new methods for the Klein-Gordon equation based on the local discontinuous Galerkin methods because LDG methods have some advantage of classical finite element methods.That is, it is easy to design LDG methods with arbitrary order of convergence, and it is suitable for mesh adaptivity and complex geometry.More importantly, DG methods can be designed to conserve mass and energy for certain systems of partial differential equations [20,21].
The outline of this paper is as follows.In Section 2, we present our semi-discrete LDG methods, prove the conservation properties of energy and linear momentum of the Klein-Gordon equation, and show the optimal convergence of our numerical solutions.All of the theoretical results from this section are based on our semi-discrete LDG schemes.In Section 3, we give rigorous proof about the superconvergence property of our numerical solution.In Section 4, we present two types of temporal discretizations and discuss their properties.In Section 5, we conclude this paper by performing some numerical experiments to test the optimal convergence, superconvergence and conservation properties.We can show that our numerical results are consistent with our theoretical results.

Semi-Discrete Local Discontinuous Galerkin Methods for the Klein-Gordon Equation
The original form of the Klein-Gordon equation is given by where φ is the wave function, c is the speed of light, m is the mass and h is the Planck constant.Applying the following transformation to Equation (1), we have the non-dimensional Klein-Gordon equation In this paper, we develop high-order numerical methods to solve Equation (3) using the method of lines.We first define the semi-discrete local discontinuous Galerkin scheme for a one-dimensional Klein-Gordon equation, and the multi-dimensional equation can be defined in the same manner.We consider the Klein-Gordon equation defined on a one-dimensional spatial domain I = [a, b] with initial conditions u(x, 0) = f (x), u t (x, 0) = g(x) and a periodic boundary condition.Before we present our numerical methods, we introduce the important notations that we will need later.We use H m to denote the L 2 -Sobolev space of order m whose equipped norm is represented by • m .We use L 2 instead of H 0 , when m = 0, and its corresponding norm is denoted by • .Similarly, we use W m,∞ to represent the L ∞ -Sobolev space of order m whose equipped norm is denoted by • m,∞ .When m = 0, we simply use • ∞ to denote the L ∞ norm.Suppose we partition the spatial domain I into N subintervals, each of which is denoted by = b be the endpoints of these subintervals, and ] for all j.For each j, let x j = 1 2 (x j− 1 2 ) be the midpoint of each subinterval I j .The mesh size is given by h := max j h j , where We define a piecewise polynomial space: V k h = {v : v| I j ∈ P k (I j ), ∀j}, where P k is the space of polynomials of degree at most k.Note that the continuity of any function in to denote the average and jump at x j+ 1

2
, respectively.In order to define local discontinuous Galerkin methods, we rewrite Equation (3) using an auxiliary variable q.Therefore, we have a (u, q) system given by u tt = q x − u and q = u x .Our local discontinuous Galerkin methods for such a (u, q) system can be defined as the following.We look for u h and q h ∈ V k h , such that depends on (q h ) + The choice of numerical fluxes plays a great role in conservation properties, stability and convergence of the numerical schemes.In this paper, we take the following numerical fluxes Alternatively, one can also choose u h = u − h , q h = q + h .

Conservation Properties
In this subsection, we show the conservation properties of our semi-discrete schemes (4) and (5).The energy and linear momentum of the Klein-Gordon equation are given by H(t) = 1 2 I u 2 t + u 2 x + u 2 dx and P(t) = I u t u x dx, respectively.In the PDE level, the energy and linear momentum of the Klein-Gordon equation are conserved.Our numerical scheme also satisfies the conservation properties.Theorem 1. (Energy Conservation) The numerical solutions u h and q h to the schemes (4) and (5) with numerical flux (6) satisfy: Proof.Let v = u h,t and q h = q − h in our semi-discrete scheme (4).We then sum over all the subintervals I j and obtain Here, we have used the periodic boundary condition of q h in the equation above.Next, we take a time derivative on both sides of Equation (5).We then let w = q h and u h = u + h and take the sum over all the subintervals I j .We can obtain Here, we have used the periodic boundary condition of u h .We take the sum of Equations ( 8) and ( 9) and get the following equation: It is easy to see that the equation above leads to the energy conservation Equation (7), and this concludes the proof.
Theorem 1 shows the energy conserving property of our numerical schemes.A direct conclusion of Theorem 1 is that u h,t 2 + q h 2 + u h 2 = C, where C is a constant that depends only on the initial data u h (x, 0), q h (x, 0) and u h,t (x, 0).Therefore, we have the following estimates: u h ≤ C, q h ≤ C and u h,t ≤ C, which implies the stability of our schemes.We then present the conservation property of linear momentum.Theorem 2. (Linear Momentum Conservation) The numerical solutions u h and q h to the schemes (4) and (5) with numerical flux (6) satisfy: Proof.Let v = q h and q h = q − h in (4) and sum over all I j , we have Let w = u h in (5) and we sum over all I j , we can get Plugging the right side of ( 13) into (12), we have We then take a t derivative on both sides of Equation ( 5), sum over all subintervals I j , and let w = u h,t , we have By taking the sum of Equations ( 14) and ( 15), we can complete the proof.
Remark 1. Theorem 2 shows that the linear momentum is conserved up to some combination of interior jumps of q h , u h and u h,t .In fact, both the linear momentum and energy can be conserved exactly, if we choose central fluxes in schemes (4) and (5), i.e., ( q h ) j+ for all j.This can be proved However, such a choice of numerical fluxes would lead to sub-optimal orders of convergence.That is, when we use piecewise-defined polynomials of degree k, the L 2 errors of u, q and u t are only of order k, instead of (k + 1) as in Theorem 3. In order to obtain the optimal accuracy, we choose alternating fluxes given by Equation (6).

Error Estimates
In this subsection, we present the error estimates of our semi-discrete schemes.Let u be the exact solution of the Klein-Gordon (3) and q = u x .We denote the error of u and q by e u := u h − u and e q := q h − q, respectively, where u h and q h are the numerical solutions to the semi-discrete schemes ( 4) and ( 5) with numerical flux (6).We then define the Gauss-Radau projections, denoted by Π − and Π + , as follows.For any function v ∈ H 1 (I), the Gauss-Radau projection of v, i.e., Π + v, is a unique function in V k h , such that for any j = 1, 2, . . ., N. Another kind of Gauss-Radau projection, Π − v, can be defined similarly.That is, for all j.Since the standard L 2 projection will also be needed, we give its definition below.The L 2 projection of v, denoted by Πv, is a unique function in V k h , such that Both Gauss-Radau and standard L 2 projections satisfy the following approximation properties.For any function v ∈ H k+1 (I), there exists a positive constant C such that where Q can be Π + , Π − or Π.Here, C depends on v k+2 when Q = Π − or Π + , and it depends on We also rewrite e q = ξ q − η q , where ξ q = q h − Π − q and η q = q − Π − q.Based on the approximation property ( 19), we know that η q and η u ≤ Ch k+1 , and are independent of the scheme.Now, we are ready to present the error estimates of our semi-discrete schemes.
Theorem 3. (Error Estimates) Let u h and q h be the numerical solutions to the schemes (4) and (5) with numerical flux (6), and u and q be the exact solutions to the Klein-Gordon Equation (3).If we take u h (x, 0) = Π + u(x, 0) and u h,t (x, 0) = Πu t (x, 0), then the following error estimates hold: where C is a positive constant that depends on t and exact solution, and is independent of the mesh size h.

Proof.
We apply numerical flux (6) to Equation ( 4) and summing over all subinterval I j , we can obtain We then take the t derivative on both sides of ( 5) and summing over all subinterval I j , we get It is easy to show that the exact solution u and q also satisfy similar equations: Subtracting ( 22) from ( 20), we have Note that we have used I η q v x dx = 0 for any v ∈ V k h , as well as the fact that ) = 0, due to the definition of Gauss-Radau projection.Similarly, subtracting (23) from (21), one can get since I η u,t w x dx = (η u,t ) + Let v = ξ u,t in (24) and w = ξ q in (25), and summing up these two equations, we have where Applying Cauchy-Schwartz inequality and the approximation property (19) to the formulation of Λ 1 , we can show that Λ 1 ≤ Ch k+1 ξ u,t + ξ q .In addition, due to the periodic boundary condition, we have Combining the estimates about Λ 1 and Λ 2 , Equation (26) leads to 1 2 Therefore, we have We integrate the equation above over time [0, t] to obtain Next, we estimate ξ u,t (•, 0) , ξ q (•, 0) and . In order to estimate ξ q (•, 0) , we consider (5) at t = 0 and get Here, we have used the fact that I q(•, 0) which leads to ξ q (•, 0) ≤ Ch k+1 .Combining all the estimates above, we have Therefore, we obtain the error estimate of u as follows: . The error estimate of q can be obtained in the same manner.

Superconvergence of Local Discontinuous Galerkin Discretization
In this section, we discuss the superconvergence property of our semi-discrete schemes.Before we present the main theorem, we give some important lemmas that will be used in the proof of the main theorem.
Lemma 1.For any x ∈ I j , let ξ u (x, t) = r j (t) + s j (x, t) x−x j h j , where r j (t) is a constant and s j (x, t) ∈ P k−1 (I j ) for any fixed t.Let s(x, t) ∈ V k−1 h , and s(x, t) = s j (x, t), ∀x ∈ I j , for j = 1, 2, . . ., N.Then, the L 2 norm of the function s(x, t) satisfies the following inequality s(•, t) ≤ Ch k+2 . (36) Proof.It is easy to show that (5) leads to Here, we have used the fact that I j η u w x dx = 0 and (η u ) + j+ 1 2 = 0, ∀j.By applying integration by parts to the term I j ξ u w x dx, we can rewrite Equation (37) as e q w dx = We then take w = s j (x, t)(x − x j+ 1 2 )/h j in (38) to obtain Combining Equations ( 39) and (40), we have We cancel the term s in the inequality above, and applying Theorem 3 to the resulting inequality, we can eventually get For convenience, let τ be any fixed time, we define τ t e u dt and E e q (x, t) = τ t e q dt.Based on these definitions, we have the following estimate.Lemma 2. For any x ∈ I j , let E ξ q = o j (t) + p j (x, t) x−x j h j , where o j (t) is a constant and p j (x, t) ∈ P k−1 (I j ) for any fixed t.Letting p(x, t) ∈ V k−1 h , and p(x, t) = p j (x, t), ∀x ∈ I j , for j = 1, 2, . . ., N.Then, the L 2 norm of the function p(x, t) satisfies the following inequality: Proof.It is easy to show that (4) leads to Here, we have used the fact that I j η q v x dx = 0 and (η q ) − j+ 1 2 = 0, ∀j.By applying integration by parts to the term I j ξ q v x dx, we can rewrite Equation (46) as We then take time integral from t to τ to obtain = 0 and .Note that the constant C in this inequality also depends on t and τ.We then cancel the term p in (52), and apply the estimates about e u,t (•, τ) , e u,t (•, t) and E e u (•, t) to the resulting inequality, we can eventually get The next lemma provides some important estimates which will be used in the proof of superconvergence.

Lemma 3. The following inequalities are satisfied
Proof.We first consider (54).Let d 3 (x) be a piecewise linear function such that d 3 (x) = (x − x j )/h j when x ∈ I j , for j = 1, 2, . . ., N. Thus, From Lemma 2, we obtain Here, we used the fact that I j r j η u,t dx = 0 in the first equality above, and applied Lemma 1 in the last inequality above.This completes the proof of (54).
Next, we consider Note that the last inequality in (57) is based on the estimates about p(•, 0) and p(•, t) from Lemma 2, as well as the approximation property about η q and E η q .Now, we are ready to present the main theorem about the superconvergence property of our scheme.
Theorem 4. (Superconvergence) Let u h and q h be the numerical solutions to the schemes (4) and (5) with numerical flux (6), and u and q be the exact solutions to the Klein-Gordon Equation (3).If we take u h (x, 0) = Π + u(x, 0) and u h,t (x, 0) = Πu t (x, 0), then the following error estimates hold: where C is a positive constant that depends on t and the exact solution, and is independent from the mesh size h.
Proof.From Equation (24), we get By rewriting I ξ u,tt v dx as d dt I ξ u,t v dx − I ξ u,t v t dx, the equation above leads to We take v = E ξ u (x, t) in the equality above and recall E ξ u (x, t) = τ t ξ u dt which implies v t = −ξ u ; thus, we can obtain From Equation ( 5), we can show We then take an integral with respect to t from τ to t on both sides of Equation ( 59), and let w = ξ q to get Note that We add Equations ( 58) and (60) and apply (61) to the resulting equation, and we can show that We take the time integral from 0 to τ in the equation above, and get the following equation Here, we have used the fact that E ξ q (x, τ) = E ξ u (x, τ) = 0, ∀x ∈ I.The first two terms on the right side of Equation ( 64) can be estimated using Lemma 3.That is, for any fixed τ, which concludes the proof of this theorem.

Time Discretizations
In Sections 2 and 3, we have defined our semi-discrete schemes and presented some conservation properties, error estimates and superconvergence properties.Here, we complete the definition of our fully-discrete schemes by providing two methods of time discretizations.That is, we consider explicit and implicit time discretizations for ( 4) and ( 5) as follows: For n = 1, 2, . .., we look for u n+1 h , q n+1 h ∈ V k h , such that Here, u n+1 h and q n+1 h represent the numerical solutions at time t n+1 := (n + 1)δt, where δt is the time step.The initialization for both schemes is given as follows: Note that Scheme 1 and Scheme 2 are explicit and implicit schemes, respectively.The property of both schemes are given in the following theorem.

Theorem 5. Let u n
h and q n h be the numerical solutions to Scheme 1 with numerical flux (6) and initialization (66) and (67); then, the following equality holds Alternatively, if u n h and q n h are the numerical solutions to Scheme 2 with numerical flux (6) and initialization (66) and (67), then the following equality holds Proof.We first prove the result for Scheme 1.We rewrite Scheme 1 as and w = q n h in the equations above, and we take the sum of the resulting equations to get Here, we have used the equality Next, we consider the result for Scheme 2. We rewrite Scheme 2 as and w = q n+1 h +q n−1 h 2 in the equations above, we can further derive It is easy to show that Equation (71) leads to which is equivalent to (69).
Remark 2. Equation (68) implies the conservation of total energy in the fully-discrete sense.Here, I u n+1 h u n h dx is an approximation of I u 2 dx, I q n+1 h q n h dx is an approximation of I q 2 dx, and approximates I u 2 t dx.

Numerical Experiments
In this section, we demonstrate the theoretical findings in Sections 2 and 3, including the optimal convergence rates of u h − u and q h − q , the superconvergence property of Π + u − u h , and conservation properties of our scheme, by some numerical tests.We consider the following initial-boundary value problem: u(x, 0) = sin(2πx), u t (x, 0) = 0, (74) The exact solution to the initial-boundary value problem is u(x, t) = sin(2πx) cos( √ 4π 2 + 1t).We apply scheme 1 defined in Section 4 to solve the problem.The numerical results by scheme 2 in Section 4 are similar, and thus we skip the discussion.
We first estimate the convergence rates of u h − u , q h − q and Π + u − u h .In these numerical experiments, we use a uniform grid with h = 1/N for various N.In order to observe the convergence order in space, we choose the time step size δt = 0.01h 2 so that the error in space dominates.We compute the numerical solutions at T = 0.5 using P k basis with k = 1, 2, 3. To estimate the convergence rate of u h − u , we take different values of h and compute u h − u .Suppose u h − u has a convergence rate of r, i.e., u h − u = Ch r where C is a positive constant; then, log 10 u h − u = r log 10 h + log 10 C. We use several points of (log 10 h, log 10 u h − u ) with different h to fit a straight line whose slope is an approximation of the convergence rate r.The convergence rates of q h − q and Π + u − u h can be obtained in the same manner.In Figure 1, we present the error of u h and q h in L 2 norm when we use P 1 basis.The empty circles represent the discrete numerical results and each red straight line represents the linear least square fitting of the empty circles.Figure 1a shows that the convergence rate of u h − u is approximately 2.0041, which is consistent with the optimal convergence rate, i.e., k + 1 with k = 1, as proved in Section 2. In addition, the results from Figure 1b verifies the convergence rate of 2 for q h − q .As for the convergence rate of Π + u − u h , it is actually greater than 2.5.Similarly, from Figures 2 and 3, we can observe the convergence rates of (k + 1) for u h − u and q h − q when P 2 and P 3 basis are used.Figure 2c shows that Π + u − u h has the convergence rate of approximately 4 when we use P 2 basis, which indicates the (k + 2) th order of convergence with k = 2. Figure 3c gives the convergence rate for the case of P 3 basis.We can see that Π + u − u h is convergent at the rate of 4.5949, which is a little over (k + 3/2) with k = 3.Based on the numerical results in Figures 1-3, as well as the discussion above, we can draw the conclusion that our numerical results coincide with our rigorous proof about error estimates and the superconvergence property.
(b) log 10 q h − q versus log 10 h.Slope = 2.1789.(b) log 10 q h − q versus log 10 h.Slope = 3.0816.Next, we check the conservation properties of our schemes using a simulation over a long time.As indicated by Theorems 1 and 2 in Section 2.2, and Theorem 5 in Section 4, the conservation properties of energy and linear momentum are independent of the degree of polynomials.Here, we only present the long time simulation when P 2 basis is used.We still consider the initial-boundary value problem (73)-(75).We choose h = 0.1, δt = 10 −4 and run the simulation up to T = 100.Using the exact solution u(x, t) = sin(2πx) cos( √ 4π 2 + 1t), we can show that the linear momentum P(t) = I u t u x dx = 0 for all t.Numerically, we approximate the linear momentum using the discrete form of P(t).That is, we compute I u n+1 h − u n h /δt • q n+1 h dx for n = 0, 1, . . ., 10 6 − 1, and present the time history of the error in linear momentum in Figure 4b.We can observe that the error oscillates from T = 0 to T = 100, and the magnitude of oscillations increase up to T = 76 and then decrease from T = 76 all the way to T = 100.Overall, the linear momentum is conserved up to the magnitude of 10 −10 .Theorem 2 implies that the rate of change for linear momentum depends on some jumps at interior interfaces, and the semi-discrete form of linear momentum is not conserved up to machine epsilon.However, our numerical results indicate that the quantity | 2 , although not equal to zero analytically, is not far away from machine epsilon numerically.Therefore, we can observe the conservation of linear momentum numerically.As for the conservation of total energy, we compute the discrete form of energy using I u n+1 h u n h dx + I q n+1 h q n h dx + u n+1 h − u n h 2 /(δt) 2 for n = 0, 1, . . ., 10 6 − 1. Figure 4a shows the error in the total energy.We can observe that the total energy is conserved up to the magnitude of 10 −9 .Therefore, we have verified the energy conservation and linear momentum conservation through the numerical experiments.

Conclusions
In this paper, high-order energy and linear momentum conserving methods are proposed for the Klein-Gordon equation.Our numerical methods are based on the local discontinuous Galerkin methods.By choosing alternating numerical fluxes (6), we can prove that the total energy of the equation is exactly conserved and the linear momentum is conserved up to some terms at element interface.Results from long time numerical simulations (up to T = 100) indicate that the linear momentum is conserved up to the magnitude of 10 −10 .Moreover, such a choice of numerical fluxes leads to optimal convergence order of u h , q h and u h,t , and the superconvergence property.It is important to mention that the schemes proposed in this paper can be directly generalized to the Klein-Gordon equation with external potential.Moreover, the results from this paper shed light on the design of high-order structure-preserving methods for nonlinear coupling systems involving the Klein-Gordon equation, for example, the Klein-Gordon-Schrödinger equations and the Klein-Gordon-Zakharov equations, which is currently under investigation.
x) be a piecewise linear polynomial on I, and d(x)| I j = x − x j+ 1 2 on each subinterval I j .It is easy to show d ∞ = h.Therefore, Equation (42) leads to log 10 Π + u − u h versus log 10 h.Slope = 3.8265.

Figure 1 .
Figure 1.log 10 E versus log 10 h, where E is the L 2 error at T = 0.5 with P 1 basis.

Figure 2 .
Figure 2. log 10 E versus log 10 h, where E is the L 2 error at T = 0.5 with P 2 basis.

Figure 3 .
Figure 3. log 10 E versus log 10 h, where E is the L 2 error at T = 0.5 with P 3 basis.

Figure 4 .
Figure 4. Time history of the error in the total energy and linear momentum with P 2 basis and h = 0.1.(left) error in the total energy; (right) error in the linear momentum.