Efficient high-order finite difference methods for nonlinear Klein – Gordon equations . I : Variants of the phi-four model and the form-I of the nonlinear Klein – Gordon equation ∗

Abstract. In this paper, the problem of solving some nonlinear Klein–Gordon equations (KGEs) is considered. Here, we derive different fourthand sixth-order explicit and implicit algorithms to solve the phi-four equation and the form-I of the nonlinear Klein–Gordon equation. Stability and consistency of the proposed schemes are studied under certain conditions. Numerical results are presented and then compared with others obtained from some methods already existing in the scientific literature to explain the efficiency of the new algorithms. It is also shown that similar schemes can be proposed to solve many classes of nonlinear KGEs.


Introduction
The nonlinear Klein-Gordon equation (KGE) where u = u(x, t) represents the wave displacement at position x and time t, and h 1 (u) is the nonlinear force, arises in the study of theoretical physics (see [18,26] and references therein).This equation is also known as the relativistic version of the Schrödinger's equation [5].It represents the equation of motion of a quantum scalar or a pseudo-scalar field (see [20] and references therein), which is a field whose quanta are spinless particles.The KGE is second order in time and therefore, it cannot be directly interpreted as a Schrödinger equation for a quantum state.Moreover, it does not admit a positive definite conserved probability density.Nevertheless, with a proper interpretation, KGE describes the quantum amplitude for finding a point particle in various places and it also describes the relativistic wave function.Moreover, complex-valued KGE is a special case of Higgs equation that is studied during interaction of scalar nucleons and mesons in particle physics [4].
Many nonlinear Klein-Gordon equations are Hamiltonian partial differential equations (PDEs), and for a wide class of force h 1 (u), they have the conserved Hamiltonian quantity (or energy) where H 1 (u) = h 1 (u).There has been enormous theoretical studies of the KGE with various potentials.There are numerous papers dealing with the existence and uniqueness of the smooth and weak solutions of these equations or studying their conservation laws (see [3,10,19,21,26,27] and references therein).In this paper, we will not study them analytically, but rather we will provide numerical procedures to obtain the solutions of the nonlinear KGE.
The solutions, thus obtained, are of interest in nuclear physics, condensed matter physics, high-energy physics, quantum field theory and many more (see [8,12,17,23] and references therein).The KGE also appears in nonlinear optics, where it plays an important role for Bose-Einstein condensates trapped in strong optical lattices formed by the interference patterns of laser beams.The KGE with vector and scalar potentials can only be solved exactly for certain potentials such as the Coloumb, Morse, or the harmonic oscillator cases, for example.In this context, it provides a good description of the scalar particles.In particular, for the case of pionic bound states in the Coloumb field of the nuclei, the KGE gives a good first order approximation of the problem [2].
In the present work, we deal with the numerical approximation of two of the best known nonlinear Klein-Gordon equations, the phi-four and its variant, the form-I of these nonlinear equations: The numerical solutions to the nonlinear Klein-Gordon equation have received considerable attention in the literature.Different low-order finite difference schemes for approximating the nonlinear KGE were discussed in [13] and [24], while the stability of the methods derived in the former paper was later analyzed in [15].Many methods have proposed spline collocation methods (see [14,17] and references therein).Radial basis functions have also been proposed for similar problems in [9].Using pseudo-spectral methods (of spectral accuracy in space and second-order accuracy in time) the solution of Klein-Gordon-Schrödinger equations were approximated in [1].The numerical solution of the damped nonlinear Klein-Gordon equations with Dirichlet boundary condition using variational method and finite element approximation was studied in [25].In [6], a numerical method based on collocation nodes is developed to solve some nonlinear Klein-Gordon equations when h 1 (u) is a polynomial.However, most of the previous methods are low-order in time or space and some of them present some limitations when they are extended to solve a large class of nonlinear Klein-Gordon equations.In this article, we will derive fourth-and sixth-order methods for one-dimensional nonlinear KGEs.The numerical methods are stable, either explicit or implicit, and easy to program.
The organization of the paper is as follows.In Section 2, we will derive higher-order (fourth-and sixth-order) explicit and implicit finite difference algorithms for the phi-four equation without forces f (x, t).We will demonstrate numerical stability when the solution of the problem is smooth enough.In Section 3, we will show that these procedures can "easily" be generalized for the form-I of the nonlinear Klein-Gordon equation.Later, in Section 4, the reader will be able to check the efficiency of the algorithms proposed in this paper, compared with others that have appeared recently in the scientific literature.In addition, we will show that these schemes can be used with small modifications to solve the PDEs previously considered with forces in space and time.

The phi-four equation
Consider the following one-dimensional nonlinear Klein-Gordon equation: with initial conditions and Dirichlet boundary conditions Let us consider a uniform mesh of the region Ω = [0, l] × [0, t e ].We can take t e = 1 and l = 1 without loss of generality.The vertices of the mesh will be (x i , t n ), where x i = ih, i = 0, 1, . . ., M , t n = nk, n = 0, 1, . . ., N , k = ∆t, and h = ∆x, and we will use the notation works exactly for the equation u tt − u xx = 0 when h = k.Using this fact and cancelling the first terms of the local truncation error, a family of high order methods was derived in [16] for the non-homogeneous hyperbolic equation u tt − u xx = q(x, t).We are using a similar idea to construct high-order finite difference methods for the phi-four equation first and later for the form-I of the nonlinear Klein-Gordon equation.http://www.mii.lt/NA

First iterations
The numerical schemes derived in this paper have three or more levels, hence it is necessary to use another procedure to obtain an approximation for the first step.Since we will derive sixth-order methods later in this section, here we will also build a sixthorder approximation for the first iteration to solve the phi-four equation ( 2) with initial conditions (3) and (4): Whenever we use a method with four or more levels, we will use the method proposed above in the first iteration and later the fourth-order three-level scheme (9) (derived below) until the other algorithm can be used (for example, if we want to use the algorithm given by equation (10), we would need to use (7) and later (9) in the second iteration).
In these cases, the order of accuracy of the algorithm is given by the following wellknown theorem.
Theorem 1. (See [22].)If the initialization of a multistep scheme uses schemes of order of accuracy r to compute the initial solution values v j for j from 1 to J such that kh r is O(h r ), and the initial data is in H ρ , where [r, ρ] is the order of accuracy of the multistep scheme, then the order of accuracy of the solution is r.

Fourth-order methods
The first method analyzed in [16] was the fourth-order scheme (for the nonhomogeneous hyperbolic equation u tt − u xx = q(x, t)) If we approximate second derivatives of q (in what follows in this paper, q = q(x, t, u) = u − u 3 ) in both space and time by three-term central methods, we obtain the formula which is a fourth-order stable algorithm for the phi-four equation as it is demonstrated in the following theorem.
Theorem 2. Formula (9) (when k = h) applied to equation ( 2) is fourth-order and unconditionally von Neumann stable when u, the solution of (2), is smooth enough.
Proof.Let us call Using Taylor series with ( 9), and after some manipulation of the formula, we obtain where Since the two first terms are zero for equation (2), the numerical scheme is fourthorder and the local truncation error of the method is 720 .
The von Neumann stability analysis is not rigorously justified for nonlinear equations, but it is often justified approximately, assuming that the solution u(x, t) (and its numerical counterpart) does not vary too rapidly, which happens with most of these nonlinear Klein-Gordon equations.When analyzing similar nonlinear equations (see, for example, the analysis in [11] for Burger's equation), the equation and its numerical method are linearized and small terms are dropped, as a small deviation between two solutions of a nonlinear partial differential equation satisfies a linearization of that PDE on the background of an exact solution.
Note that −144 + 144 cos(θ) 2 0, but still, for some values of α and θ, the radical can be positive.However, the radical can be considered as )h (we consider h 1 since we are proving stability) and |g 1 ||g 2 | = 1.Thus, in any case, both roots satisfy |g i | 1 + C(u 0 )h and the numerical method is von Neumann stable.
In [7], it is established that this condition is not sufficient for stability.In addition, it is required that the powers of the matrix where β = −2(4αh 2 + (12 + αh 2 ) cos(θ))/(−12 + αh 2 ), are determined by its eigenvalues.A simple calculus give us Remark 1.It is not surprising that all the roots of the amplification polynomial of method (9) satisfy the condition In fact, as we stated before, all the numerical schemes in this paper come from modifying (8), which is well-known to be stable (with both roots of the amplification polynomial |g i | = 1), adding to the scheme terms that are O(k 2 ) uniformly in θ.That is, since |g i | 1 when k = 0, considering in addition that the solution of the PDE is smooth, if we use Taylor series, we can conclude that for small values of k, the roots of the amplification polynomials of the schemes in this paper satisfy |g i | 1 + C(u 0 )k (see Corollary 2.2.2 and Theorem 2.2.3 in [22]), and only two eigenvalues are ≈ 1.
This suggests that when the solution u of ( 2) is smooth enough, it is possible to use any of the following options to demonstrate stability: (i) a similar procedure as in the previous Theorem, or (ii) a similar procedure as in Theorems 2.2.1 and 2.2.3, and Corollary 2.2.2 in [22].Hence, in what follows whenever we analyze a method, we will only study the local truncation error.
The first numerical method, presented in (9), was obtained approximating the second derivatives of q in both space and time by three-term central methods.Thus, ( 9) is a threelevel algorithm, where it is necessary to solve a third-order equation at each point for the phi-four model.Although this is relatively easy to calculate, for other nonlinear Klein-Gordon equations, this could be a little more difficult, even though a simple Newton-Raphson scalar method can be used with good results.
In case, we want to use a fully explicit scheme, we would only need to approximate the term (q tt ) n i in ( 8) by the three-term backward scheme )/h 2 (and (q xx ) n i by the central scheme in space as for ( 9)) to obtain a third-order four-level method.If we want to obtain a fourth-order explicit algorithm, then it is necessary to use five levels, which we can achieve by replacing (q tt ) n i in ( 8) by h 2 (and again replacing (q xx ) n i by the central scheme in space).Theorem 3. The four-level explicit formula (when k = h), applied to equation (2), is only third order and its local truncation error is the following: .
Theorem 4. The explicit formula (when k = h), applied to equation (2), is fourth order and its local truncation error is the following: .
http://www.mii.lt/NAAnother possibility is considering q(x, t, u) = u − u 3 .In that case, a simple calculus provide us the formula Here, we can consider central schemes in space (for u x and u xx ) and both central and backward schemes in time (for u t ).In this way, we obtain three different algorithms similar to those derived above: one explicit with order three, a fourth-order three-level implicit scheme, and a fourth-order five-level explicit method.

Sixth-order methods
In a similar way as in the previous subsection, we can derive sixth-order methods from the sixth-order algorithm derived in [16] for the wave equation.However, in this case, it is necessary to approximate the term (q xx ) n i + (q tt ) n i with a fourth-order accurate formula, and the term (q xxxx ) n i + (q ttxx ) n i + (q tttt ) n i with a second-order method.Thus, we can obtain, for example, the following two sixth-order algorithms.
Theorem 5.The local truncation error of the sixth-order six-level implicit method (when k = h), applied to equation (2), can be written briefly as Theorem 6.The local truncation error of the sixth-order explicit and seven-level scheme (when k = h), applied to equation (2), can be written briefly as Obviously, methods given by equations ( 12) and ( 13) can be used only when i = 2, . . ., M − 2, so it is necessary to calculate other schemes for the cases i = 1 and i = M − 1.For i = 1, it is necessary to approximate (q xx ) n i and (q xxxx ) n i with analogous formulae but evaluated in q n 0 , q n 1 , q n 2 , q n 3 , q n 4 , and q n 5 .We need to proceed similarly when i = M − 1, approximating (q xx ) n i and (q xxxx ) n i with the formulae evaluated in q n M , q n M −1 , q n M −2 , q n M −3 , q n M −4 , and q n M −5 .It is also possible to construct similar sixth-order methods using the fact that q(x, t, u) = u − u 3 , calculating in terms of u and its derivatives in space and time, and finally replacing these derivatives with fourth-order numerical derivative formulae.
With both procedures it is also possible to obtain higher-order methods for the phifour equation.

Nonlinear Klein-Gordon equation (form-I)
As we commented above, the phi-four equation is only a particular case of the form-I of the nonlinear Klein-Gordon equation again with initial conditions (3) and ( 4), and boundary conditions ( 5) and (6) as above.Hence, procedures as those employed in the previous Section can be used again, and actually the numerical schemes are very similar to those obtained in the previous theorems.
Thus, for the first iteration, now the approximation is Fourth-order numerical schemes, equivalent to those given by equations ( 9) and ( 11), and whose properties were studied in Theorems 2 and 4, are given with the following theorem.
Theorem 7. The implicit three-level formula and the five-level explicit method are both fourth order, and their local truncation error are respectively Numerical stability of the algorithm given by ( 16) can be proven as in (2), just changing u − u 3 by u − u m and writing α = 1 − u m−1 0 instead of α = 1 − u 2 0 .Again, it is necessary that the solution of problem ( 14) is smooth enough.For the stability of the numerical scheme (17) it is better to proceed as it was suggested in Remark 1.
Analogously, it is also relatively easy to construct sixth-order and eight-order methods, and other higher-order schemes with a similar procedure as above.Actually, Theorems 5 and 6 can be used for the problem given by formula (14), it is only necessary to write q = u − u m instead of q = u − u 3 .

Numerical examples
In this section, we offer some numerical examples developed with the methods shown in the preceding sections.

First numerical example
For the first numerical example, we shall consider problem (1) with the traditional initial and boundary conditions ( 3)-( 6) such that the solution is ) .
As we commented above for our fourth-and sixth-order methods, we need to approximate the solution at the first step, and for that we use equation (7) to calculate it.In this numerical experiment, we will employ the implicit fourth-order method (9) and the explicit sixth-order method (13).For the former one, we require additional approximations of the solution at the first levels.Hence, there we can employ (9) at the first steps.
As we can check in Table 2, this methods converge at the expected rate.There, we show the values r i = log ci err i /err 1 , where err i is the error obtained with h i and As we also would like to study the numerical convergence of the sixth-order method, and for this one, we would need to employ several times another algorithm at the first steps, we will proceed in two different ways to better study the convergence rate of the Table 1.Errors for the fourth-order method (9), 1st example.scheme: (i) we will use the real solution for the first steps and (ii) we will employ the implicit fourth-order method (9) in these first steps.First, in Tables 3 and 4 the numerical errors and the study for the convergence rate are shown with the real solution in the first steps.Later, in Tables 5 and 6, we proceed in a similar way with procedure (ii).The convergence rate is as expected (see Theorem 1).
As we can check in Tables 4 and 6 this method converges at the expected rate.There, we showed the values r i = log ci err i /err 1 , where err i is the error obtained with h i and c i = h i /h 1 .

Second numerical example
For the second numerical example we shall consider a problem with unknown solution, where the equation and initial conditions are as follows: and the following boundary conditions: As in the previous example for our fourth-and sixth-order methods, we need to approximate the solution at the first step, and for that we use equation (15) to calculate it.In this numerical experiment we will employ the implicit fourth-order method ( 16) for iterations 2 and 3, and the explicit fourth-order scheme (17) in the followings steps.
We shall also consider a sixth-order method, in this case, the one given by formula (13).As iterations 2 to 5 are unknown, we have used ( 16) before we can employ the sixth-order explicit formula (13).
As in the previous problem, we can check in Table 8 that this method converges at the expected rate.Again, we showed the values r i = log ci err i /err 1 , where err i is the error obtained with h i and c i = h i /h 1 .
We also would like to study the numerical convergence of the sixth-order method and for this one, we need to employ several times a fourth-order method at the first steps as we comment before.In Tables 9 and 10, the numerical errors and the study for the convergence rate are shown.

Third numerical example
In this test example, we compare our schemes with the fourth example used by Rashidinia et al. in [17].That paper is devoted to the numerical solution of a one-dimensional nonlinear Klein-Gordon equation, which is given in the following form: where s(t, x) = −2 + x 2 cosh(t + x) + x 6 cosh 3 (t + x) − 4x sinh(t + x), and subjected to the initial conditions with the exact solution u(x, t) = x 2 cosh(x + t).They obtained the boundary conditions from the exact solution, and they solved this problem with the algorithm derived in [17] with h = 0.01 and k = 0.0001.
As the reader can check, this equation has not been exactly considered previously in the paper (since the phi-four equation is u tt − u xx − u + u 3 = 0, a sign is different and s(t, x) = 0 in this case).Also the interval, where the equations are considered, is different, ranging from −1 to 1 instead of 0 to 1.However, similar algorithms can be derived and again stability for (smooth enough) functions can be proved as in previous methods from Remark 1 and Duhamel's principle (a scheme is stable for the equation P k,h v = f if it is stable for the equation P k,h v = 0 as it is shown in [22,Chap. 9]).
In this case, we have used a scheme similar to the one given by equation (7) for the first iteration, and later, for the following 4 iterations, we have used an implicit method similar to the one derived in equation (8).Finally, we have applied a six-order explicit algorithm similar to the one given by equation (13).We have used the same problem than they have, but in this our with h = k = 0.01.Hence, our k is much bigger than theirs, which means a large reduction in computing time.However, the methods studied here are higher-order.We did not want to compare our schemes with theirs in other numerical problems since it looks "unfair", but we considered that we had to do it in one instance to make clearer the differences to the reader.Table 11 offers the comparison results.We calculated the errors at t e = 1 with the infinity and L 2 norm and also the root mean square (RMS) which can be calculated in this case as the error with the L 2 norm divided by the root square of the number of nodes (201).We took the values with the other methods from Table 4 in [17].

Conclusions
In this paper, the problem of solving some nonlinear Klein-Gordon equations (KGE) has been considered.More exactly, we have derived different fourth-and sixth-order explicit and implicit algorithms to solve the phi-four equation and the form-I of the nonlinear Klein-Gordon equation.We have also demonstrated that similar schemes can be obtained for similar equations.The stability and consistency of the proposed schemes were studied under certain conditions of smoothness of the solution of these equations.
We have provided numerical results to explain the efficiency of the new algorithms, their stability, and the convergence rate.We have compared the results obtained with the new methods against others obtained with some schemes already existing in the scientific literature.In comparison, they show an improvement in both, the computational effort and the error rate, as they are smaller than in the other schemes.
Regarding the future work, we think it would be interesting to derive and study the numerical implementation of numerical schemes for other nonlinear Klein-Gordon when G(u) in u tt − u xx − G(u, u x , u t , . ..) = 0 is not a polynomial or suffers perturbations.