THIRD-ORDER EXPONENTIAL INTEGRATOR FOR LINEAR KLEIN–GORDON EQUATIONS WITH TIME AND SPACE-DEPENDENT MASS

. Allowing for space-and time-dependence of mass in Klein–Gordon equations resolves the problem of negative probability density and of violation of Lorenz covariance of interaction in quantum mechanics. Moreover it extends their applicability to the domain of quantum cosmology, where the variation in mass may be accompanied by high oscillations. In this paper we propose a third-order exponential integrator, where the main idea lies in embedding the oscillations triggered by the possibly highly oscillatory component intrinsically into the numerical discretisation. While typically high oscillation requires appropriately small time steps, an application of Filon methods allows implementation with large time steps even in the presence of very high oscillation. This greatly improves the efficiency of the time-stepping algorithm. Proof of the convergence and its rate are nontrivial and require alternative representation of the equation under consideration. We derive careful bounds on the growth of global error in time discretisation and prove that, contrary to standard intuition, the error of time integration does not grow once the frequency of oscillations increases. Several numerical simulations are presented to confirm the theoretical investigations and the robustness of the method in all oscillatory regimes


Introduction
In this paper we consider the Klein-Gordon equation  2  2 (, ) = ∆(, ) + (, )(, ),  ∈ [ 0 ,  ],  ∈ T  (,  0 ) =  0 (),   (,  0 ) =  0 (), with periodic boundary conditions 1 .Here (, ) denotes the unknown wave function that we wish to approximate numerically and −(, ) > 0 stands for the square of the time-and space-dependent mass, which may be presented or approximated by the truncated, modulated Fourier expansion.Thus we assume in this paper applications and is consistent e.g. with finite difference, finite element and spectral approaches.In our numerical examples we assumed that the problem, consistently with (1.1), is equipped with periodic initial and boundary conditions and used a Fourier spectral collocation in space.
Outline of the paper.In Section 2 we recall basic properties of the Duhamel formula and indicate how they can be adapted to Klein-Gordon-type equations.Section 3 presents the comprehensive analysis of the convergence of the method.In Section 4 we illustrate our theoretical findings with numerical experiments.In addition, we compare various existing methods with the new approach and highlight the favorable behavior of the new method with respect to high frequencies   .

Derivation of the method
As already stated in the introduction, we are concerned solely with time discretisation: in practical application it should be of course accompanied by a space discretisation, but our framework is consistent with an arbitrary choice of the latter.For this reason in the further part of this manuscript we suppress dependance on space variable  and write (),   (), () instead of (, ),   (, ), (, ).In order to derive third order methods we assume that  ′′′ ∈  1 ([ 0 ,  ]).In this setting we can express (1.1) as the following abstract evolution equation with formal solution given by Duhamel formula, where , and  •  = −∆.
Note that, −∆ being positive definite, we can also choose a positive definite .
In computational implementation we resort to full discretisation and then functions (),   (), () are replaced by vectors and the unbounded operator  by a suitable positive2 , finite-dimensional operator   .
Let us assume that  /ℎ =  and   = ℎ.Iterating Duhamel's formula we obtain the following relations at the time steps  +1 =   + ℎ,  = 0, . . .,  − 1: As in the case of Gautschi-type methods [9], the above representation of the solution is the foundation of our numerical scheme.
Remark 2.1.In this manuscript we understand that  ′ (  ) and  ′′ (  ) are the first and second, respectively, time derivatives of function ( • , ) at point   .In particular we write  ′ ( 0 ) instead of  0 .
The system (2.3) reduces to two time-stepping formulas for the values of the functions () and  ′ () at the time step  +1 , namely ) Table 1.Algorithm Ξ [3] (time integration) for finding the approximate solution   and  ′  in the time interval [ 0 ,  ] at points   =  0 + ℎ, ℎ = ( −  0 )/ .Important part of this pseudocode are the integrals.In the case when they cannot be computed analytically, we recommend Filon-type methods described in the further part of this paper.Practical approximation of the integrals is presented in Appendix A. Discretisation in space is general and depends on boundary conditions.In the numerical examples, presented in Section 4, we assume periodic initial and boundary conditions and use Fourier collocation for space discretisation.
Numerical scheme Ξ [3] ]︁  end do The main difficulty of the above formulation is that it is an implicit system of equations.Given the initial conditions we resort within the integral to the Taylor expansion where Remark 2.2.The above approximation could be directly incorporated into (2.4) and (2.5) because  ′′ (  ) can be immediately recovered from the initial condition (for  = 0) and from the original problem (1.1) for  ≥ 1.The latter, however would lead to additional numerical evaluations.For this reason for  ≥ 1 we resort again to the Taylor formula where Ē = Remark 2.3.We further observe that the errors   ( ) and Ē contain the third derivative of the unknown function , which obviously depends on oscillations   .Indeed, due to the structure of the equation (1.1) expression  3   involves   , so we formally have that Given that  3   grows to infinity together with ω, it is natural to expect that   ( ) and Ē would also grow to infinity and that these values would affect the overall error constant of approximation.However, as will be shown in Lemma 3.4 the error of the method does not depend on the oscillations.Application of (2.6) and (2.7) in (2.4), (2.5) leads directly to a numerical scheme Ξ [3] presented in Table 1, where   and  ′  are numerical approximations of (  ) and (  ), respectively.
Remark 2.4.If the integrals appearing in numerical scheme Ξ [3] cannot be computed analytically, they need to be approximated.It is important to emphasize that the high performance of numerical approach Ξ [3] and its third-order convergence are obtained only once we approximate the possibly highly-oscillatory integrals (containing (  +  )) with a quadrature of order at least 4 (locally) where the error constant does not depend on the oscillatory parameters   and where no relation between ℎ and ω is required.For this reason standard quadratures like Gauss-Legendre are unequal to the task (because their error constants depend on the time derivatives of the integrants, in particular on the highly oscillatory part) and we need a more specialised approach.Here we propose to use the Filon-type methods which obtain desired order in time step ℎ, whose error constants do not depend on (possibly) large ω and where no relation between ℎ and ω need be imposed.
Filon methods are based on the one-dimensional oscillatory integral, where   () is an -oscillatory function, by which we broadly mean a function that oscillates rapidly in .Its organising principle is to approximate the non-oscillatory function  () by a polynomial () such that For Filon-type methods we refer the reader to [8] and [11], where they have been introduced and comprehensively analysed.In our case it is enough to use the quadratic polynomial, Remark 2.5.In the scheme Ξ [3] we apply a (locally) fourth-order Filon-type method to the integrals of the form ∫︀ ℎ 0  ()e  .It is easy to verify that for () defined as in (2.8) we have where  is independent of  ≫ 1.This, in turn, means that where the integral on the right hand side can be evaluated explicitly by integration by parts.
Note that, unless a Filon method (or similar highly oscillatory quadrature) is used, the method will fail unless the time step ℎ > 0 is minute, essentially ℎ < ω−1 .This phenomena is illustrated in Example 4.1, see Figure 1.

The convergence of the method
It is not straightforward to conclude convergence and stability from the equations in scheme Ξ [3] , because naive calculations of global error lead to exponential growth in number of time steps taken in the interval [ 0 ,  ].For this reason we need to reformulate the equations (2.4) and (2.5)., where results have been obtained with Ξ [3] method with applications of Filon-type method and high order Gauss-Legendre quadratures.
Proof.We proceed by induction.One can easily check that the formulas for ( 1 ) and  ′ ( 1 ) obtained by (3.1) and (3.2) are equivalent to these obtained by (2.4) and (2.5), respectively.Let us assume that the same holds for (  ) and  ′ (  ) for certain 1 being positive definite, we can apply the trigonometric identities, and conclude with This proves equivalence of the values of function  obtained by (2.4) and (3.1).A similar result for  ′ obtained by (2.5) and (3.2) can be shown in a similar manner.
′′′ (),  = 0, . . .,  − 1, Ē0 = 0, and The global errors   := (  ) −   and  ′  :=  ′ (  ) −  ′  ,  = 1, . . ., , of method Ξ [3] can be expressed with formulas where Proof.Based on the Theorem 3.1 and Remark 3.2 the global errors   and  ′  can be expressed as Recalling the initial condition, ( 0 ) =  0 ,  ′ ( 0 ) =  ′ 0 , and the fact that  ′′ 0 is restored from the equation (cf. the first line of scheme Ξ [3] ), we observe that Next, recalling that Ē0 = 0 and letting for  = 1, . . .,  and  = 0, . . .,  − 1, we easily conclude that Given that sin(( and subtracting consecutive global errors (3.7) and (3.8) their growth can be represented in the form where, telescoping we derive the claim of the lemma.Proof.Functions ℛ  and ℛ ′  ,  = 1, . . ., , depend on the errors   ( ) and Ē , which in turn depend on the third derivative of the unknown function, The second term of the above expression involves multiplications by   , which require careful attention.For this reason we write   ( ) =  0  ( ) +    ( ),  = 0, . . ., , where and Ē0 = 0, Ē = Ē0  + Ē  ,  = 1, . . ., , where It is trivial to observe, that  0  ( ) ∼ We represent Ē  in a similar manner.This leads to the conclusion that    ( ) ∼  Next we focus our attention on the structure of ℛ  and ℛ ′  .Obviously Moreover we observe that ‖2 −1 sin )︀ ‖ ≤ ℎ‖‖.Equipped with this estimates we easily conclude the claim of the Lemma.Remark 3.5.Obviously the unboundedness of operator  may raise concerns, because of the presence of ‖‖ in the estimate.We recall, though, that the theme of this paper is time integration and that eventually space discretization must be combined with our narrative.In practice, this means that in applications  will be approximated with some bounded operator   , in particular that ‖  ‖ will be bounded.
Finally we are ready to formulate and prove the theorem on the convergence of our method.
Proof.Let  := max ∈[0, ] |()|.Based on the Theorem 3.3 and Lemma 3.4 we obtain for  ≥ 2 the following inequalities Let now   and  ′  be the upper bounds of ‖  ‖ and ‖ ′  ‖ respectively,  = 0, . . ., , therefore and for  ≥ 2 We conclude that and that the following inequalities are satisfied: Let us recall now that ‖‖ ≥ 1.Based on (3.11)-(3.13)we observe that ℛ ≤ ℛ,   ≤  ′  which in turn means that and where  depends on  , , max This means that global error of the method ‖( ) −   ‖ is bounded uniformly in ω by at least ℎ 2 and that in the limit as ℎ → 0 we obtain third order global error.

Numerical examples
We thereby focus on two equations.In Example 4.1 we consider the wave equation with a highly oscillating term of single frequency  and discuss behaviour of numerical methods with regard to the size of this parameter.We are concerned with the effectiveness of scheme Ξ [3] depending on the type of quadratures used for approximation of the integrals appearing in the scheme.Moreover we compare various numerical approaches and show that each of them is effective in a different oscillatory regime.In Example 4.2 we consider the wave equation with multiple frequencies and run comparisons of different time stepping methods.It will be seen that the proposed third-order method performs exceedingly well in this setting.
In all the experiments we assume that the solution is confined to the spatial interval [ 0 ,   ] and assume periodic boundary conditions.We divide the interval into  = 200 sub-intervals of length ∆ = (  −  0 )/ and discretise spatial operators with Fourier collocation method, described in detail in [14].As a reference solution we apply to the semi-discretised equations a 6th order method based on self-adjoint basis of Munthe-Kaas & Owren [3] with a step size ℎ = 10 −6 .We carry out numerical integration for presented method using various time steps and measure the  2 absolute error at the final time  .The errors are plotted in doublelogarithmic scale.In Figure 1 parameter  takes values 103 (on the left), 10 4 (in the centre) and 10 5 (on the right), respectively.The solution to the equation (4.1) is approximated with the Ξ [3] method, with the integrals computed with the fourth-order Filon method (as presented in Rem.2.5).The derived solution is compared against the results of the Ξ [3] method complemented by Gauss-Legendre quadratures of fourth, sixth and eighth order, respectively.As can be seen, even high orders of Gauss-Legendre quadratures fall well short of the third (globally) order method unless the time step ℎ is significantly smaller than the frequency .This confirms our claim that Filon methods (or other highly oscillatory quadrature methods) need be applied in the presence of high oscillation.
In Figure 2 we compare the method Ξ [3] against other schemes from the literature, like well known 2nd and 4th order Runge-Kutta methods (RK [2] and RK [4] ).We will also consider the 4th-order method Σ [4] 3 from [1] (denoted here as BBCK [4] ) which is designed for non-oscillatory potentials, and the 3rd-order asymptotic method from [6] (denoted by Asympt [3] ), appropriate for extremely high oscillations 3 .The comparisons are analysed in three regimes, with  = 10,  = 500 and  = 1500.As can be seen method Asympt [3] proves it efficiency only in presence of high oscillations and fails in the case of  = 10.All other methods work very well in slowly oscillatory regime.Second order Runge-Kutta method loses its effectiveness already for  = 500.It is clear from the middle and right columns that the second and fourth order methods RK [2] ,RK [4]  and BBCK [4] require time step ℎ to be smaller than  −1 .The new method Ξ [3] maintains third global error and proves its effectiveness in all (from slowly to highly) oscillatory regimes and does not require the time step ℎ to be smaller than  −1 .We can also observe, that the error constant does not grow with the growth of .According to our computations the error of the new method Ξ [3] scales like ℎ 3 globally uniformly in , see Figure 3.Note that the error constant is not affected by ω = 10 5 , a very high frequency, and that the new method behaves well for all time steps.Moreover, as expected, the asymptotic method performs very purely.This is caused by the small frequency  1 = 1 which appears in the input term.On the other hand the large oscillations like  6 = 10 5 sabotage methods that approximate poorly highly oscillatory integrals.
Our two examples demonstrate vividly the value of our proposed method Ξ [3] in comparison with both classical and asymptotic methods.

Figure 1 .
Figure 1.Accuracy and time of computation in seconds of the numerical solution to problem (4.1),where results have been obtained with Ξ[3] method with applications of Filon-type method and high order Gauss-Legendre quadratures.

Figure 2 .
Figure 2. Comparison of accuracy of methods known in the literature, where three oscillatory regimes in equation 4.1 are considered.