Numerical schemes for general Klein–Gordon equations with Dirichlet and nonlocal boundary conditions

. In this work, we address the problem of solving nonlinear general Klein–Gordon equations (nlKGEs). Different fourth- and sixth-order, stable explicit and implicit, ﬁnite difference schemes are derived. These new methods can be considered to approximate all type of Klein– Gordon equations (KGEs) including phi-four, forms I, II, and III, sine-Gordon, Liouville, damped Klein–Gordon equations, and many others. These KGEs have a great importance in engineering and theoretical physics. The higher-order methods proposed in this study allow a reduction in the number of nodes, which might also be very interesting when solving multi-dimensional KGEs. We have studied the stability and consistency of the proposed schemes when considering certain smoothness conditions of the solutions. Additionally, both the typical Dirichlet and some nonlocal integral boundary conditions have been studied. Finally, some numerical results are provided to support the theoretical aspects previously considered.


Introduction
We have studied in this paper the general Klein-Gordon equations, i.e., where q(t, x, u) is a sufficiently differentiable function and t e ∈ R + . Without loss of generalization, l = 1 can be considered.
c Vilnius University, 2018 Typical initial conditions are utilized and Dirichlet and/or nonlocal integral conditions are studied u(1, t) = l 0 ψ(x, t)u(x, t) dx + g 2 (t), 0 t t e .
We will consider that q, f ,f , g 1 , g 2 , φ, and ψ are known functions, and we will assume enough smoothness in these functions to get the desire orders of convergency of the algorithms. Furthermore, it is assumed that those functions satisfy those conditions so that the solution to equations (1)- (5) exists and is unique. The KGEs are a family of partial differential equations (PDEs) that arises when solving relativistic quantum mechanics and field theory problems. These equations model many different phenomena in engineering (including acoustics or biomedical engineering as it can be seen, for example, in [7,17,19,29]), theoretical physics in general, and high energy physics in particular [18].
Many nlKGEs are Hamiltonian partial differential equations, and for a wide class of force h 1 (u), they have the following conserved Hamiltonian quantity (or energy): where ∂Q(x, t, u)/∂u = q(x, t, u). We can find in the scientific literature an extensive research on the KGE with various potentials. Most of these studies are related to the existence and uniqueness of the smooth and weak solutions of these equations and sometimes to their conservation laws (see [6,13,32,33,39,41] and references therein). However, this paper does not deal with the analytical solution, but with the numerical procedures that are necessary in order to obtain efficient solutions to the nlKGE.
Although numerical solutions for nlKGEs have received considerable attention in the literature, most of these schemes are low-order. Different finite difference schemes for approximating nlKGE were discussed in [21] and [36]. Stability of the methods derived in the former paper was later analyzed in [24]. Some articles have proposed spline collocation methods (see [22,30] and references therein). Radial basis functions have also been proposed for similar problems in [12]. Klein-Gordon-Schrödinger equations were approximated in [3] using pseudo-spectral methods (with spectral accuracy in space and second-order accuracy in time). Numerical solutions of the damped nonlinear Klein-Gordon equations with Dirichlet boundary conditions using a variational method and a finite element approximation were studied in [37]. In [7], a numerical method based on collocation nodes was developed to solve some Klein-Gordon equations when q(x, t, u) is a polynomial.
In [14,15], we developed some higher-order schemes for some Klein-Gordon equations, phi-four and forms I, II, and III of the nlKGEs. In this paper, we show how to extrapolate these techniques to all types of nlKGEs, including variants and perturbed forms I, II, and III, but also sine-Gordon, double sine-Gordon, sinh-Gordon, Liouville, etc.
However, none of these articles study these PDEs together with nonlocal integral conditions. Many mathematical models arising in various fields of science and engineering can be expressed as partial differential equations with nonlocal initial or boundary conditions. For example, PDE-based models are common in different physical processes and also in thermoelasticity, ecology, chemistry, semiconductor modelling, medicine, and biotechnology [1,[8][9][10][11] when it is impossible to determine the boundary or initial values of the unknown function. In addition to that, they are also used in problems related to control or inverse PDEs (see [23,38]). Nonclassical boundary and initial-boundary value problems with integral and/or discrete nonlocal boundary conditions were widely studied for numerous equations. In the last decades, numerical methods for the solution of partial differential equations with nonlocal conditions have been developed and analyzed very actively (e.g., see [2, 4, 5, 25-27, 31, 35, 40] and references therein). In this work, we have analyzed equation (1) with standard initial condition (typical Dirichlet), but also nonlocal integral conditions (4)-(5) in boundary.
Thus, we have developed the numerical schemes in Section 2: Subsection 2.1 presents the first step of the methods, while Subsection 2.2 describes the process to derive highorder finite difference schemes for nlKGEs and their convergence properties, and Subsection 2.3 explains how nonlocal boundary conditions can be solved. Later, in Section 3, numerical results associated to different types of problems are shown to demonstrate the efficiency of the new schemes, but also to numerically support the theoretical results previously given. Finally, some conclusions are given in Section 4.

Numerical schemes
The numerical schemes developed in this paper consist on three or more levels, therefore, it is necessary to use another procedure to obtain an adequate approximation for the first step. Later in this section, we will derive the explicit and implicit fourth-and sixth-order finite difference schemes before addressing the boundary conditions.
As it is usually, in the case when constructing finite difference schemes, we have considered a uniform mesh of the region Ω = [0, l] × [0, t e ]. The vertices of the mesh will be (x i , t n ), where t n = nk, n = 0, 1, . . . , N , with k = ∆t and x i = ih, i = 0, 1, . . . , M = l/h, with h = ∆x. As in [14,15], we will require that k = h to obtain higher-order schemes.

First step
For the first iteration, we merely consider to build a sixth-order approximation. u 0 i and (u t ) 0 i are known from the initial conditions (2) and (3). Additionally, since

becomes a different equation depending on the nonlinear Klein-Gordon equation and the values of
, and so on.
, and so on. For some full examples on these first iterations, readers can obtain full expressions for the phi-four equation, form-I, form-II, and form-III of the nonlinear KGE in [14,15].

Fourth-and sixth-order methods
Once we have the first iteration, we can employ the fourth-and sixth-order methods. In order to obtain them, we start from the algorithm which exactly solves the equation u tt − u xx = 0. Now, let us introduce the term q(t, x, u) in the right-hand side of the former equation, and let us operate to obtain the Taylor series of the new equation. The two first terms of the truncation error obtained in this way can be approximated through three-term central finite difference methods (i.e., the traditional schemes for second derivatives of q). Thus, the fourth-order implicit scheme (for the nonhomogeneous hyperbolic equation (1) with h=k) is derived in the following way: It is possible to check that the leading term of the truncation error is However, if we use backward differences in time to approximate the second derivative of q, then we can develop the following fourth-order explicit scheme whose leading term of the truncation error is It is possible to obtain sixth-order methods for the general nlKGEs in a similar way as it is described in [27]. We can start from the formula and approximate the term (q xx ) n i + (q tt ) n i with a forth-order accurate formula and the term (q xxxx ) n i + (q ttxx ) n i + (q tttt ) n i with a second-order method. Thus, it is possible to obtain the following two sixth-order algorithms: [27].) The local truncation error of the sixth-order six-level implicit method (when k = h) can be written briefly as Theorem 2. (See [27].) The local truncation error of the sixth-order explicit and sevenlevel scheme (when k = h) can be written briefly as It is clear that the methods given by equations (13) and (14) 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 the first case, 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 also need to proceed similarly when i = M − 1 approximating (q xx ) n i and (q xxxx ) n i with the symmetrical 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 .

Stability of the numerical schemes
So far, we have developed several methods with fourth-and sixth-order local truncation errors. In order to demonstrate convergence of the proposed algorithms, we now need to show that they are also stable. In this subsection, we have collected some theoretical results explained in [34], and we will use them to prove the stability of the numerical algorithms derived in previous sections. An important application of Fourier analysis is the von Neumann analysis of stability of finite difference schemes. Let us assume that the Fourier inversion formula is given by and the Fourier transform to the solution of the scheme in the following time step is equivalent to multiplying the Fourier transform of the solution by the so called amplification factor g(hξ, k, h) (which is normally a function depending on θ = ξh, h, and k): u n+1 (ξ) = g(hξ, k, h)û n (ξ).
Therefore, we have the following results: These theorems are devoted to the first-order hyperbolic equation when the employed methods are one-step algorithms. When the equation is second-order and the algorithms are not one-step, we can consider the following result:  Φ(g, θ) for a secondorder time-dependent equation is explicitly independent of h and k, then the necessary and sufficient condition for the finite difference scheme to be stable is that all roots, g ν (θ), satisfy the following conditions: (i) |g ν (θ)| 1; (ii) if |g ν (θ)| = 1, then g ν (θ) must be at most a double root. Proof. The idea for this proof is analogous to the one employed in Theorem 7 in [14] for some numerical schemes proposed for some specific nonlinear Klein-Gordon equations.
To the best of our knowledge, the von Neumann stability analysis has not been 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 nonlinear Klein-Gordon equations. When analyzing similar nonlinear equations (see, for example, the analysis in [16] for Burger's equation), the equation and its numerical method are linearized, and small terms are dropped (when k, h are small enough) 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. In what follows, we are assuming that k is small enough, so the equation and its numerical method can be linearized for the stability analysis.
We now remind the readers that all the numerical schemes proposed above, (8), (10), (13), and (14), are modifications of (7) consisting only in the addition of terms that are O(k 2 ). Therefore, the amplification factors of the new schemes are modifications result in the addition of terms that are O(k 2 ) uniformly in ξ of the amplification factor of the original method (7).
Therefore, an analogous proof to Theorem 4 can be used. But, in this case, the equation is second-order in time. Hence, it is necessary that modifications are O(k 2 ) additions to the original algorithm. Thus, instead of satisfying Theorem 3, the original method has to satisfy Theorem 5.
Since the amplification polynomial Φ(g, θ) associated to (7) is their roots, g 1,2 (θ) = e ±iθ , satisfy |g ν (θ)| = 1. Hence, both conditions are satisfied. https://www.mii.vu.lt/NA The stability for PDEs with nonlocal boundary conditions is much more complex in general (see [20,31] and references therein). As we will check in the next section, when φ(x, t) = 0 and/or ψ(x, t) = 0, then it is necessary to approximate u n+1 0 and/or u n+1 M through equations (15) and/or (16), and therefore, the difference scheme depend on φ(x, t) and ψ(x, t). Actually, we consider that stability for these Klein-Gordon equations with nonlocal conditions cannot be studied with the previous amplification polynomials and should be studied through spectral structure of some transition matrices, which complicate much more this procedure.
Unfortunately, we feel that it is impossible to study, with only one theorem, the stability for all kind of nonlinear Klein-Gordon equations with all the methods proposed in this paper and for all nonlocal boundary conditions.
If one wants to study the stability for one particular Klein-Gordon equation (fixed q(x, t, u)) and for a given numerical scheme (any of the methods proposed in this paper), then we would propose, first, to fix φ(x, t) and ψ(x, t). However, there are still many factors in the transition matrices affected by the k-value (k = h also appears in (15) and (16) multiplying φ(x, t) and ψ(x, t) as it is shown below). Hence, we think that it is necessary to take several fixed values of k = h and numerically study the stability for these particular values and test.

Nonlocal boundary conditions
In [28], a parabolic problem with nonlocal boundary conditions was studied. Although in this work we are solving a nonlinear hyperbolic equation, we can use the same procedure to obtain the linear equations for u 0 and u M .
In case we want to derive fourth-order approximations, we use fourth-order Simpson's composite formula, which allows us to obtain and where and In case we want to consider the sixth-order schemes, we can also develop them as it is explained in [28]: and where and 3 Numerical examples In this section, we will compare the new schemes proposed in the previous sections to show the properties studied above.
To the best of our knowledge, nonlocal boundary conditions together with nonlinear Klein-Gordon equations are numerically analyzed for the first time in a scientific paper.
We computed the numerical solutions for ε = 0 and ε = 0.3 with the four schemes. When we considered ε = 0, the four methods provided errors ∼ O(10 −14 ) − O(10 −16 ) (round-off errors). Thus, we were able to check that these codes are exact for many examples where solutions are low-order polynomials.
Logically, they are also very accurate when solutions are approximately some of these polynomials. In Fig. 1(c), errors for the four schemes are compared when ε = 0.3.

Conclusions
In this work, we have developed and analyzed fourth-and sixth-order stable explicit and implicit finite difference schemes for general nonlinear Klein-Gordon equations. Stability and consistency of these schemes were studied theoretically and numerically for the case of nonlocal boundary conditions together with nonlinear Klein-Gordon equations. We have included for this purpose three numerical examples in which we clearly appreciate the different convergence rate and stability of these schemes. It was shown that they can be considered for a wide range of different nlKGEs, including with typical Dirichlet, but also some nonlocal integral boundary conditions.