An Energy Conserving Numerical Scheme for the Klein–Gordon Equation with Cubic Nonlinearity

: In this paper, we develop a numerical scheme that conserves the discrete energy for solving the Klein-Gordon equation with cubic nonlinearity. We prove theoretically that our scheme conserves not just discrete energy, but also other energy-like discrete quantities. In addition, we prove the convergence and the stability of the scheme. Finally, we present some numerical simulations to demonstrate the performance of our energy-conserving scheme.


Introduction
The Klein-Gordon equation is a well-studied equation in mathematical physics in terms of radiation theory, nonlinear optics and general relativity of scattering [1][2][3].In particular, much work has been carried out on the Klein-Gordon equation with respect to wave collisions and resonance behavior [4].Recently, we considered the Klein-Gordon equation with quintic nonlinearity from the analytical point of view and were able to obtain a number of new wave solutions [5].In this paper, our focus is on the numerical study of the initial value problem of the Klein-Gordon equation with cubic nonlinearity with the initial conditions where α and β are arbitrary positive constants, and u 0 (x) and u 1 (x) are known smooth functions.Here, one should note that the non-local version of (1), where the dependence of previous time history is considered, gives a time-fractional nonlinear Klein-Gordon equation [6].There have been a number of analytical and numerical studies to understand the solutions of nonlinear time and space-fractional Klein-Gordon equations [7][8][9][10][11].The numerical studies have focused on various techniques to discretize the fractional derivatives.However, no special attention was given to the nonlinearities in the equations, as we do in Section 2.
In order to construct a numerical scheme for solving (1), we re-cast the initial value problem as the following initial boundary value problem: with the initial conditions u(x, 0) = u 0 (x), u t (x, 0) = u 1 (x), and the boundary conditions u(x L , t) = 0, u(x R , t) = 0, 0 ≤ t ≤ T.
Here x L is negative, and x R is positive with |x L | and x R being large values so that the finite spatial domain x L ≤ x ≤ x R mimics the infinite domain −∞ < x < ∞.Moreover, T can be large.In some cases, the boundary conditions (4) may be replaced by The boundary conditions at u(x L , t) and u(x R , t) for the problem (2) correspond to the asymptotic conditions for u(x, t) of (1) as x goes to −∞ and ∞.The solution u(x, t) of the initial boundary value problem (2)-( 4) formally satisfies the following energy identity: In the literature, one can find a number of numerical schemes with conservation properties for solving the nonlinear Klein-Gordon equation.For example, a three-level finite difference method that conserves energy was developed in [12] , while other finite-difference algorithms that preserve energy or linear momentum were studied in [13].In addition, there are schemes that were constructed using a variational iteration method [14,15], a homotopyperturbation idea [16,17], radial basis functions [18], spline-collocation approach [19] and discrete Fourier transforms [20] to solve the Klein-Gordon equation under various conditions.Further, some of the recent studies on the nonlinear Klein-Gordon equation have involved making use of pseudo-spectral discretization methods [21], employing a differential quadrature method with cubic B-splines [22] and domain decomposition methods [23].Other studies include making use of the sinc-collocation idea along with a discrete gradient method to study the Klein-Gordon-Schrödinger equation [24] and developing a higher order method for the Klein-Gordon equation employing a local discontinuous Galerkin method [25].However, in [25], the numerical simulations were carried out only for the linear Klein-Gordon equation.In [26], energy-preserving schemes were constructed for higher dimensional Klein-Gordon equations using the discrete gradient method and Duhamel principle.In addition, there have been a couple of interesting analytical studies of the Klein-Gordon equation, one utilizing an operational matrix method with clique polynomials [27] and the other a series method using differential transforms [28].
Most of the existing numerical methods investigate the conservation of discrete energy only numerically.If one is to validate the numerical results of an energy conserving numerical scheme, it is important to prove theoretically that the scheme conserves the discrete energy.The work in [29] carried out a theoretical study of four explicit finite difference schemes for solving the Klein-Gordon equation.In the spirit of [29], this paper presents an implicit conservative finite difference scheme for the initial boundary value problem (2)-( 4).It should be noted that in [30] implicit finite difference schemes were studied for the coupled system of Klein-Gordon-Zakharov equations.Later, more work along the same lines was conducted in [31] for the same system.Even though there is some similarity, in contrast to those works, our study considers not just the conservation of the discrete energy, but other energy-like discrete quantities as well.A predictor-corrector idea is employed to deal with the nonlinearity which appears in the problem.Furthermore, we give some a priori estimates and then prove by the discrete energy method that the difference scheme is stable and second-order convergent.Some numerical results are presented to illustrate the theoretical results.Three-dimensional plots are displayed to demonstrate the sensitivity of the discrete energy and other discrete quantities to the choices of time steps, wave speed, and coefficients α and β.

Finite Difference Scheme and Its Conservation Law
Before we propose the conservative finite difference scheme for the Klein-Gordon equation with cubic nonlinearity (2)-( 4), we give some notations as follows when we discretize the space and time domains: where h and τ are the step sizes of space and time, respectively.In addition, we define the following inner product and norms: It should be noted that in the following, C stands for a general positive constant that may take different values on different occasions.In addition, for brevity, we omit the subscript 2 of w n 2 .
Lemma 1.For any two mesh functions {w m } and {v m }, m = 0, 1, 2, . . ., M, there is the identity This lemma can be easily proved using the notational definitions directly.Let U n m be the difference approximation of u(x, t) at (x m , t n ); that is, U n m ≈ u(x m , t n ).In addition, assume that u 0 (x m ) = U 0 (x m ) and u 1 (x m ) = U 1 (x m ).Now, we consider the following finite difference scheme for the Klein-Gordon equation with cubic nonlinearity (2)-( 4): In order to employ the finite difference scheme (6), we need initial values at two different time levels.They are chosen from the initial conditions given in (3) such that making use of a fictitious time level −1.
The boundary conditions are as below.
It should be pointed out that our two-time-level split approximation of the nonlinear cubic term is very different than the standard nonlinear approximation.As will be seen later in Theorems 1 and 2, this split approximation makes the theoretical analysis easier.
In (6), the explicit forms of (U n m ) t t and (U n m ) x x are given as follows: As noted before, the solution u(x, t) of the initial boundary-value problem (2)-(4) satisfies the following energy identity: Now, we present some properties of our finite difference scheme.
Theorem 1.The difference scheme (6)-( 8) possesses the following property: where Proof.Computing the inner product of (6) with U n − U n−1 , we have In the computation of the above equation, we have used the boundary conditions and Lemma 1. Now, using the Taylor's series expansions for u(x m , t n+1 ) and u(x m , t n−1 ) about u(x m , t n ), we can easily show that Here u n m = u(x m , t n ) and ∂ ∂t is denoted by .So, if the higher order terms of τ are neglected beyond τ 4 , we have Therefore, for the finite difference approximation U n m of u n m , we obtain the relationship In a similar fashion, we can obtain Using these relationships for the finite difference approximations, we obtain where This completes the proof of the theorem.
Theorem 2. The difference scheme (6)-( 8) possesses the following invariant: where Proof.Computing the inner product of (6) with U n+1 − U n−1 , we have In the computation of Equation ( 13), we have used the boundary conditions and Lemma 1.By adding and subtracting 1  2 U n x 2 , α 2 U n 2 , and 4 to the left-hand side of Equation (13) and rearranging the terms, we obtain Hence, result (11) follows from Equation (14).This completes the proof.Now, from Equations (10) and (12), we can easily observe that and therefore, it follows from (9) that Moreover, from Equation (15) and Equation (11), we have or equivalently, In addition, from Equation (9), we obtain Therefore, it follows from (17) and (18) that Theorem 3. The difference scheme (6)-( 8) possesses the following property: where Moreover, if E n is given by Equation (12), we have Proof.From (16) and (19), we have In addition, from Equations (22), and (11), we obtain Hence, result (20) follows from Equation (24), and result (22) follows from Equation (23).This completes the proof.
It should be pointed out that even though our numerical scheme (6) (with two-timelevel split approximation) is second order, it does not immediately follow that every discrete quantity that will be conserved will also be conserved up to second order.As we have shown in Theorems 1-3, if a discrete quantity, such as Q n , U n t 2 , or Ẽn is defined using only the n th time level, then each one of them will be conserved up to order one.On the other hand, the discrete energy E n defined at two time levels n and (n + 1) is shown to be conserved (Theorem 2) without any order restrictions.

Some a Priori Estimates for the Numerical Solutions
In this section, we will obtain some a priori estimates for the numerical solutions of the scheme (6).Our work makes use of the lemmas presented in [32].
Lemma 3 (Gronwall's Inequality).Suppose that the nonnegative mesh functions {w(n where B l (l = 1, 2, . . ., N) are nonnegative constants.Then, for any 0 ≤ n ≤ N, there is then the following estimates hold: Proof.From Equation (11), we have In addition, from Equation (25), we obtain Therefore, it follows from the last two equations that Hence, we obtain In addition, we can obtain the following estimate by Lemma 2: This completes the proof.

Convergence and Stability of the Difference Scheme
In this section, we will discuss the convergence and the stability of the difference scheme (6)- (8).First, we define the truncation error by Lemma 4. Assume that the conditions of Theorem 4 are satisfied, and u(x, t) ∈ C 4,4 , then the truncation error of the difference scheme (6)-( 8) satisfies By Taylor's expansion, Lemma 4 can be proved directly.Moreover, we note that the approximation of the initial condition (7) has the truncation error of order O(τ 2 ), which is consistent with the scheme.Now, we are going to analyze the convergence of the difference scheme (6)- (8).Let us set Theorem 5. Assume that the conditions of Lemma 4 are satisfied.Then the solution of the difference scheme (6)-( 8) converges to the solution of the problem stated in (2)-( 4) with order O(τ 2 + h 2 ) in the L ∞ norm for U n .
Proof.Subtracting (6) from (26), we obtain Then computing the inner product of (27) with e n+1 − e n−1 , we have where Using Young's inequality ab ≤ 1 4 a 2 + b 2 , we have Therefore, where  then by (28) and Lemma 4, we have

So, we have
Summing (29) up for n (1 ≤ n ≤ N), we obtain and hence, we have Applying Gronwall's inequality (Lemma 3), we obtain Therefore, we obtain From the discrete initial conditions, we know that e 0 and e 1 are of second-order accuracy, then Hence, the following inequalities can be obtained by (30): This completes the proof of Theorem 5.
It should be remarked that since our boundary value problem (2)-( 4) involves second derivatives of u(x, t) in time and space, in order for the difference scheme (6)-( 8) to be a consistent second order method in both time and space, foundational theory in numerical analysis dictates that u(x, t) ∈ C 4,4 (also, see [31]).If for example, u(x, t) ∈ C 3,3 , still the finite difference scheme (6)-( 8) works, but now, it will be a consistent first order method in both time and space, i.e., u(x, t) ∈ C 4,4 is not an essential condition for the method to be consistent.
In the same way as above and under the conditions of Theorem 5, we can also prove that the solution U n of the difference scheme (6)-( 8) is stable in the sense of norm .∞ .

Numerical Results
In this section, we will test the efficiency of our numerical scheme by considering a number of simulations.A predictor-corrector idea is employed to deal with the nonlinearity which appears in the problem.
Let us first define the "error" functions as

Bell Solitary Wave Solution
We consider the following initial boundary value problem of the Klein-Gordon equation with cubic nonlinearity subject to the initial conditions and the homogenous Dirichlet boundary conditions Note that the initial conditions are derived from the exact solitary wave solution of (31) given by [33] This exact solution is of bell shape and represents a soliton which travels with velocity c and whose amplitude is 2α β .The exact solution u(x, t) of the above initial boundary value problem satisfies the following energy identity: This is fairly straightforward and is obtained by applying Equation (32) in Equation (5) at t = 0.For our computations, we consider parameters α = 1, β = 1 π 2 , and c = 1.5.Hence, the approximate value of the constant E is Since (6)-( 8) is a three-time-level numerical scheme, in order to get the computer simulation started, at the beginning, we need initial values at two different time levels.For our computations of bell solitary wave solution and kink solitary solution, respectively, these initial values were obtained from u(x, 0) and u t (x, 0) making use of the respective exact solutions given by (32) and (33).
In Figure 1, the solitary wave computed by the numerical scheme (6)-( 8) is compared with the wave of exact solution at time T = 5.As one can see, both waves are indistinguishable-the numerical solution simply overlaps the exact solution.The curves of discrete energy E n and discrete quantities Ẽn , Q n , and U n t 2 obtained by the numerical scheme (6)-( 8) at a larger T value (T = 7) are plotted in Figure 2. The figure shows that the numerical scheme (6)-( 8) possesses very good conservation properties when compared to the theoretical results.In order investigate the influence of the time-step size τ, the computations were repeated with a fixed space step h = 0.1 and a different time-step size τ = 0.0125.Figure 3 shows the sensitivities of E n , Ẽn , Q n , and U n t 2 to the time-step size τ.We can easily see that the discrete quantities Ẽn , Q n , and U n t 2 are more sensitive than the discrete energy E n to the changing of the time-step size τ.
Figure 4 shows the comparison between the exact solution and the numerical solution with h = 0.1 and τ = 0.0125 at time T = 50 for c = 0.1.One can easily see that the solitary wave solution computed by the numerical scheme (6)-( 8) agrees very well with the exact solution.In addition, the curves of discrete energy E n and discrete quantities Ẽn , Q n , and U n t 2 obtained by the numerical scheme (6)-( 8) are plotted in Figure 5.This shows that the numerical scheme (6)- (8) possesses extremely good conservation properties.Table 1 gives the numerical errors for the scheme (6)- (8) with different h and τ at time T = 50.In fact, the errors are presented for mesh widths h and time steps τ as they are halved.Using simple arithmetic, one can easily verify that the L ∞ error decreases as second order in time and space when τ and h are halved.Tables 2 and 3 show the conservation of discrete energy E n and discrete quantities Ẽn , Q n , and U n t 2 computed by the numerical scheme (6)- (8) with h = τ = 0.1 at time T = 10, 20, 30, 40, and 50.Moreover, Table 4 gives the errors between exact and approximate discrete energies and quantities with different velocities at different times in the case when h = 0.1 and τ = 0.0125.∞ and e n with c = 0.2, h = 0.1, and τ = 0.05, 0.025, 0.015, and 0.01 at time T = 1.The error functions are computed at different values for α and β.Hence, for a small velocity c, the number of error oscillations decreases as α decreases and β increases-i.e., when the cubic term dominates the linear term in the Klein-Gordon nonlinearity.∞ and e n computed when α equals 0.8, 0.9, 1.0, 1.1, and 1.2 and β equals 0.8, 0.9, 1.0, 1.1, and 1.2.

Conclusions
In this paper, we constructed a finite difference scheme that conserves the discrete energy (and some other discrete quantities) for solving the Klein-Gordon equation with cubic nonlinearity.Theoretical analysis is provided to show the conservation properties of the numerical scheme.In addition, we obtain theoretical error estimates and prove the stability and the convergence of the scheme.Finally, we carry out a number of computer simulations using the scheme.In particular, we consider the cases where the solutions are either traveling pulses or traveling wave fronts.The numerical simulations demonstrate that our method performs very well in both instances-conserving the discrete energies and producing accurate and stable solutions.One observation is that if it is imperative to conserve the other discrete quantities along with the discrete energy, one may have to choose a smaller time step.This is because since the conservation of the discrete quantities are correct up to the order of the spatial mesh and the time step, at instances, some of the discrete quantities, other than the discrete energy, are susceptible to an increasing time step.However, this does not affect the performance of our method.One can still conserve the discrete energy and obtain excellent numerical results that are stable and accurate.In addition, in the case of traveling wave fronts with low speeds, we find that our scheme performs well (with no error oscillations) if the cubic term is dominant compared to the linear term (i.e., larger β and smaller α).As we noted in the introduction, there are a few energy conserving explicit finite difference schemes in the literature for solving the Klein-Gordon equation.However, because of the explicitness, the stability of these schemes is conditional resulting in restrictive choices for the spatial mesh width and time step.In contrast, since our energy-conserving scheme is an implicit scheme, the stability is unconditional, and we do not have any restrictions on the spatial mesh width or the time

Figures 6 and 7
Figures 6 and 7 show the error functions e n∞ and e n with c = 0.2, h = 0.1, and τ = 0.05, 0.025, 0.015, and 0.01 at time T = 1.The error functions are computed at different values for α and β.Hence, for a small velocity c, the number of error oscillations decreases as α decreases and β increases-i.e., when the cubic term dominates the linear term in the Klein-Gordon nonlinearity.

Table 2 .
Discrete energy E n and discrete quantities Ẽn , Q n , and U n t 2 with h = τ = c = 0.1.

Table 3 .
Discrete energy E n and discrete quantities Ẽn , Q n , and U n t 2 with h = τ = 0.1 and c = 0.2.

Table 4 .
The errors between exact and approximate discrete energies and quantities with different velocities at different times.