Long time error analysis of finite difference time domain methods for the nonlinear Klein-Gordon equation with weak nonlinearity
aa r X i v : . [ m a t h . NA ] M a r Long time error analysis of finite difference time do-main methods for the nonlinear Klein-Gordon equa-tion with weak nonlinearity
Weizhu Bao , Yue Feng ∗ , Wenfan Yi Department of Mathematics, National University of Singapore, Singapore 119076 School of Mathematics and Econometrics, Hunan University, Changsha, 410082,Hunan Province, P. R. China
Abstract.
We establish error bounds of the finite difference time domain (FDTD) meth-ods for the long time dynamics of the nonlinear Klein-Gordon equation (NKGE) witha cubic nonlinearity, while the nonlinearity strength is characterized by ε with 0 < ε ≤
1a dimensionless parameter. When 0 < ε ≪
1, it is in the weak nonlinearity regime andthe problem is equivalent to the NKGE with small initial data, while the amplitudeof the initial data (and the solution) is at O ( ε ) . Four different FDTD methods areadapted to discretize the problem and rigorous error bounds of the FDTD methodsare established for the long time dynamics, i.e. error bounds are valid up to the timeat O ( ε β ) with 0 ≤ β ≤
2, by using the energy method and the techniques of eitherthe cut-off of the nonlinearity or the mathematical induction to bound the numericalapproximate solutions. In the error bounds, we pay particular attention to how errorbounds depend explicitly on the mesh size h and time step τ as well as the small pa-rameter ε ∈ ( ] , especially in the weak nonlinearity regime when 0 < ε ≪
1. Our errorbounds indicate that, in order to get “correct” numerical solutions up to the time at O ( ε β ) , the ε -scalability (or meshing strategy) of the FDTD methods should be takenas: h = O ( ε β /2 ) and τ = O ( ε β /2 ) . As a by-product, our results can indicate error boundsand ε -scalability of the FDTD methods for the discretization of an oscillatory NKGEwhich is obtained from the case of weak nonlinearity by a rescaling in time, while itssolution propagates waves with wavelength at O ( ) in space and O ( ε β ) in time. Ex-tensive numerical results are reported to confirm our error bounds and to demonstratethat they are sharp. AMS subject classifications : 35L70, 65M06, 65M12, 65M15, 81-08
Key words : nonlinear Klein-Gordon equation, finite difference time domain methods, long timeerror analysis, weak nonlinearity, oscillatory nonlinear Klein-Gordon equation.
DedicatedtoProfessorJieShenontheoccasionofhis60thbirthday ∗ Corresponding author.
Email addresses: [email protected] (W. Bao), [email protected] (Y. Feng), [email protected]
Consider the nonlinear Klein-Gordon equation (NKGE) with a cubic nonlinearity on atorus T d ( d = ) [23, 27, 36, 37] as ∂ tt u ( x , t ) − ∆ u ( x , t )+ u ( x , t )+ ε u ( x , t ) = x ∈ T d , t > u ( x ,0 ) = φ ( x ) , ∂ t u ( x ,0 ) = γ ( x ) , x ∈ T d . (1.1)Here t is time, x ∈ R d is the spatial coordinates, u : = u ( x , t ) is a real-valued scalar field, 0 < ε ≤ φ ( x ) and γ ( x ) are two given real-valued functionswhich are independent of ε . The NKGE is a relativistic (and nonlinear) version of theSchr ¨odinger equation and it is widely used in quantum electrodynamics, particle and/orplasma physics to describe the dynamics of a spinless particle in some extra potential[4, 7, 13, 22, 33, 34, 36]. Provided that u ( · , t ) ∈ H ( T d ) and ∂ t u ( · , t ) ∈ L ( T d ) , the NKGE (1.1)is time symmetric or time reversible and conserves the energy [5, 19], i.e., E ( t ) : = Z T d (cid:20) | ∂ t u ( x , t ) | + |∇ u ( x , t ) | + | u ( x , t ) | + ε | u ( x , t ) | (cid:21) d x ≡ Z T d (cid:20) | γ ( x ) | + |∇ φ ( x ) | + | φ ( x ) | + ε | φ ( x ) | (cid:21) d x : = E ( ) = O ( ) , t ≥
0. (1.2)We remark here that, when 0 < ε ≪
1, rescaling the amplitude of the wave function u by introducing w ( x , t ) = ε u ( x , t ) , then the NKGE (1.1) with weak nonlinearity can bereformulated as the following NKGE with small initial data, while the amplitude of theinitial data (and the solution) is at O ( ε ) : ∂ tt w ( x , t ) − ∆ w ( x , t )+ w ( x , t )+ w ( x , t ) = x ∈ T d , t > w ( x ,0 ) = εφ ( x ) , ∂ t w ( x ,0 ) = εγ ( x ) , x ∈ T d . (1.3)Again, the above NKGE (1.3) is time symmetric or time reversible and conserves theenergy [5, 19], i.e.,¯ E ( t ) : = Z T d (cid:20) | ∂ t w ( x , t ) | + |∇ w ( x , t ) | + | w ( x , t ) | + | w ( x , t ) | (cid:21) d x = ε E ( t ) ≡ Z T d (cid:20) ε | γ ( x ) | + ε |∇ φ ( x ) | + ε | φ ( x ) | + ε | φ ( x ) | (cid:21) d x : = ¯ E ( ) = O ( ε ) . (1.4)In other words, the NKGE with weak nonlinearity and O ( ) initial data, i.e. (1.1), isequivalent to it with small initial data and O ( ) nonlinearity, i.e. (1.3). In the following,we only present numerical methods and their error bounds for the NKGE with weaknonlinearity. Extensions of the numerical methods and their error bounds to the NKGEwith small initial data are straightforward. There are extensive analytical results in the literature for the NKGE (1.1) (or (1.3)).For the existence of global classical solutions and almost periodic solutions as well asasymptotic behavior of solutions, we refer to [10–12,15,40–42] and references therein. Forthe Cauchy problem with small initial data (or weak nonlinearity), the global existenceand asymptotic behavior of solutions were studied in different space dimensions andwith different nonlinear terms [25, 26, 31, 35, 38]. Recently, more attentions have beendevoted to analyzing the life-span of the solutions of the NKGE (1.3) [25, 32]. The resultsindicate that the life-span of a smooth solution to the NKGE (1.3) (or (1.1)) is at least upto time at O ( ε − ) [16, 18]. For more details related to this topic, we refer to [17, 21] andreferences therein.For the numerical aspects of the NKGE (1.1) (or (1.3)), different numerical methodshave been proposed and analyzed in the literatures [5, 14, 20, 44], including the finite dif-ference time domain (FDTD) methods [5, 14, 20, 44], exponential wave integrator Fourierpseudospectral (EWI-FP) method [5, 6, 9], multiscale time integrator Fourier pseudospec-tral (MTI-FP) method [4], etc. In these results, the error bounds are normally valid up tothe time at O ( ) . Since the life-span of the solution of the NKGE (1.1) can be up to thetime at O ( ε − ) , it is a natural question to ask how the performance of a numerical methodfor (1.1) up to the time at O ( ε − ) , i.e. long time error analysis. In other words, one has toestablish error bounds of the numerical method for (1.1) up to the time at O ( ε − ) insteadof the classical error bounds which are only valid up to the time at O ( ) . The purpose ofthis paper is to carry out rigorous error analysis of four widely used FDTD methods forthe NKGE (1.1) in the long time regime. In our error bounds, we pay particular attentionto how the error bounds depend explicitly on the mesh size h and time step τ as well asthe small parameter ε ∈ ( ] . In our numerical analysis, besides the standard techniqueof the energy method and the inverse inequality, we adapt the cut-off of the nonlinearityfor the conservative methods, and resp., the mathematical induction for nonconservativemethods, to obtain a priori bound of the numerical solution in the l ∞ norm. Based onour rigorous error bounds, in order to obtain “correct” numerical approximations of theNKGE (1.1) (or (1.3)) up to the long time at ( ε − β ) with 0 ≤ β ≤ ε -scalability (or meshing strategy requirement) of the FDTD methods when 0 < ε ≪ h = O ( ε β /2 ) and τ = O ( ε β /2 ) .As a by-product, by rescaling the time as t → t / ε β with 0 ≤ β ≤ O ( ) in space and O ( ε β ) in time. The FDTD methods to (1.1)and their error bounds over long time can be extended straightforwardly to the oscilla-tory NKGE up to the time at O ( ) . With the error bounds, the ε -scalability (or meshingstrategy) of the FDTD methods for the oscillatory NKGE can be drawn.The rest of the paper is organized as follows. In Section 2, different explicit/semi-implicit/implicit and conservative/nonconservative FDTD discretizations are presentedfor the NKGE (1.1) and their properties of the stability, conservation and solvability areanalyzed. In Section 3, we establish rigorous error estimates of the FDTD methods for the NKGE (1.1) over long time dynamics. Extensive numerical results are reported in Section4 to confirm our error bounds. In Section 5, we extend the FDTD methods and theirerror bounds to an oscillatory NKGE. Finally, some conclusions are drawn in Section 6.Throughout this paper, we adopt the notation p . q to represent that there exists a genericconstant C >
0, which is independent of the mesh size h and time step τ as well as ε suchthat | p | ≤ Cq . In this section, we adapt four different FDTD methods to discretize the NKGE (1.1) andanalyze their properties, such as stability, energy conservation and solvability. For sim-plicity of notations, we shall only present the numerical methods and their analysis forthe NKGE (1.1) in one space dimension (1D). Thanks to tensor grids, generalizations tohigher dimensions are straightforward and results remain valid with minor modifica-tions. In 1D, consider the following NKGE ∂ tt u ( x , t ) − ∂ xx u ( x , t )+ u ( x , t )+ ε u ( x , t ) = x ∈ Ω = ( a , b ) , t > u ( x ,0 ) = φ ( x ) , ∂ t u ( x ,0 ) = γ ( x ) , x ∈ Ω = [ a , b ] , (2.1)with periodic boundary conditions. Choose the temporal step size τ : = ∆ t > h : = ∆ x >
0, and denote M = ( b − a ) / h being a positive integer and the grid points and time steps as: x j : = a + jh , j = M ; t n : = n τ , n = X M = { u = ( u , u ,..., u M ) T | u j ∈ R , j = M , u = u M } and we always use u − = u M − and u M + = u if they are involved. The standard discrete l , semi- H and l ∞ norms and inner product in X M are defined as k u k l = h M − ∑ j = | u j | , k δ + x u k l = h M − ∑ j = | δ + x u j | , k u k l ∞ = max ≤ j ≤ M − | u j | , ( u , v ) = h M − ∑ j = u j v j ,with δ + x u ∈ X M defined as δ + x u j = ( u j + − u j ) / h for j = M − u nj be the numerical approximation of u ( x j , t n ) for j = M , n ≥ t = t n as u n = ( u n , u n ,..., u nM ) T ∈ X M . We introduce the finitedifference operators as δ + t u nj = u n + j − u nj τ , δ − t u nj = u nj − u n − j τ , δ t u nj = u n + j − u nj + u n − j τ , δ + x u nj = u nj + − u nj h , δ − x u nj = u nj − u nj − h , δ x u nj = u nj + − u nj + u nj − h .Here we consider four frequently used FDTD methods to discretize the NKGE (2.1):I. The Crank-Nicolson finite difference (CNFD) method δ t u nj − δ x (cid:16) u n + j + u n − j (cid:17) + (cid:16) u n + j + u n − j (cid:17) + ε G (cid:16) u n + j , u n − j (cid:17) = n ≥
1; (2.3)II. A semi-implicit energy conservative finite difference (SIFD1) method δ t u nj − δ x u nj + (cid:16) u n + j + u n − j (cid:17) + ε G (cid:16) u n + j , u n − j (cid:17) = n ≥
1; (2.4)III. Another semi-implicit finite difference (SIFD2) method δ t u nj − δ x (cid:16) u n + j + u n − j (cid:17) + (cid:16) u n + j + u n − j (cid:17) + ε (cid:16) u nj (cid:17) = n ≥
1; (2.5)IV. The leap-frog finite difference (LFFD) method δ t u nj − δ x u nj + u nj + ε (cid:16) u nj (cid:17) = j = M − n ≥
1. (2.6)Here, G ( v , w ) = F ( v ) − F ( w ) v − w , ∀ v , w ∈ R , F ( v ) = Z v s ds = v v ∈ R . (2.7)The initial and boundary conditions in (2.1) are discretized as u n + = u n + M , u n + − = u n + M − , n ≥ u j = φ ( x j ) , j = M , (2.8)where the initial velocity γ ( x ) is employed to update the first step u by the Taylor ex-pansion and the NKGE (2.1) as u j = φ ( x j )+ τγ ( x j )+ τ h δ x φ ( x j ) − φ ( x j ) − ε (cid:0) φ ( x j ) (cid:1) i , j = M . (2.9)It is easy to check that the above FDTD methods are all time symmetric or time re-versible, i.e. they are unchanged if interchanging n + ↔ n − τ ↔ − τ . In addition,the LFFD (2.6) is explicit and might be the simplest and the most efficient discretizationfor the NKGE (2.1) with the computational cost per time step at O ( M ) . The others areimplicit schemes. Nevertheless, the CNFD (2.3) and SIFD1 (2.4) can be solved via eithera direct solver or an iterative solver with the computational cost per time step depend-ing on the solver, which is usually larger than O ( M ) , especially in two dimensions (2D)and three dimensions (3D). Meanwhile, the solution of the SIFD2 (2.5) can be explicitlyupdated in the Fourier space with O ( M ln M ) computational cost per time step, and suchapproach is valid in higher dimensions. Let T > ≤ β ≤
2, and denote σ max : = max ≤ n ≤ T ε − β / τ k u n k l ∞ . (2.10)Following the von Neumann linear stability analysis of the classical FDTD methods forthe NKGE in the nonrelativistic limit regime [5, 29], we can conclude the linear stabilityof the above FDTD methods for the NKGE (2.1) in the following lemma. Lemma 2.1. (linear stability) For the above FDTD methods applied to the NKGE (2.1) up to thetime t = T ε − β , we have:(i) The CNFD (2.3) is unconditionally stable for any h > τ > and < ε ≤ .(ii) When h ≥ , the SIFD1 (2.4) is unconditionally stable for any h > and τ > ; and when < h < , this scheme is conditionally stable under the stability condition < τ < h / p − h , h >
0, 0 < ε ≤
1. (2.11) (iii) When σ max ≤ ε − , the SIFD2 (2.5) is unconditionally stable for any h > and τ > ; andwhen σ max > ε − , this scheme is conditionally stable under the stability condition < τ < p ε σ max − h >
0, 0 < ε ≤
1. (2.12) (iv) The LFFD (2.6) is conditionally stable under the stability condition < τ < h / q + h ( + ε σ max ) , h >
0, 0 < ε ≤
1. (2.13)
Remark 2.1.
The stability of schemes (2.5) - (2.6) is related to σ max , dependent on theboundedness of the l ∞ norm of the numerical solution u n at the previous time step. Theconvergence estimates up to the previous time step could ensure such a bound in the l ∞ norm, by making use of the inverse inequality, and such an error estimate could berecovered at the next time step, as given by the Theorems presented in Section 3.For the CNFD (2.3) and SIFD1 (2.4), we can show that they conserve the energy in thediscretized level with the proofs proceeding in the analogous lines as those in [5, 30, 37]and we omit the details here for brevity. Lemma 2.2. (energy conservation) For n ≥ , the CNFD (2.3) conserves the discrete energy asE n : = k δ + t u n k l + n + ∑ k = n k δ + x u k k l + n + ∑ k = n k u k k l + ε h M − ∑ j = h ( u nj ) +( u n + j ) i ≡ E . (2.14) Similarly, the SIFD1 (2.4) conserves the discrete energy as ˜ E n : = k δ + t u n k l + h M − ∑ j = ( δ + x u nj )( δ + x u n + j )+ n + ∑ k = n k u k k l + ε h M − ∑ j = h ( u nj ) +( u n + j ) i ≡ ˜ E , n ≥
0. (2.15)
Based on Lemma 2.2, we can show the unique solvability of the CNFD (2.3) at eachtime step as follows.
Lemma 2.3. (solvability of CNFD) For any given u n , u n − (n ≥ ), the solution u n + of theCNFD (2.3) is unique at each time step.Proof. Firstly, we prove the existence of the solution for the CNFD (2.3). To simplify thenotations, we denote the grid function [[ u ]] n ∈ X M with [[ u ]] nj = u n + j + u n − j j = M , n ≥
1. (2.16)For any u n − , u n , u n + ∈ X M , we rewrite the CNFD (2.3) as [[ u ]] n = u n + τ F n ([[ u ]] n ) , n ≥
1, (2.17)where F n : X M → X M with F nj ( v ) = δ x v j − (cid:20) + ε ( | u n − j | + | v j − u n − j | ) (cid:21) v j , j = M , n ≥
1. (2.18)Define a map K n : X M → X M as K n ( v ) = v − u n − τ F n ( v ) , v ∈ X M , n ≥
1. (2.19)It is obvious that K n ( n ≥
1) is continuous from X M to X M . Moreover, the fact ( K n ( v ) , v ) = k v k l − ( u n , v )+ τ (cid:20) k δ + x v k l + k v k l + ε (cid:16) | u n − | + | v − u n − | , v (cid:17)(cid:21) ≥ ( k v k l −k u n k l ) k v k l , n ≥
1, (2.20)implies lim k v k l → ∞ ( K n ( v ) , v ) k v k l = ∞ , n ≥
1. (2.21)Then, we can conclude that there exists a solution v ∗ such that K n ( v ∗ )= k u n k l + k δ + x u n k l ≤ E n = E , n ≥
0. (2.22)Hence, by employing the discrete Sobolev inequality [2, 39], we can obtain k u n k l ∞ . k u n k l + k δ + x u n k l . √ E , n ≥
0. (2.23)
For any v ∈ X M , we define a functional S ( v ) : X M → R as S ( v ) : = M − ∑ j = " − u nj + u n − j τ − δ x u n − j + u n − j + ε (cid:16) u n − j (cid:17) v j + M − ∑ j = (cid:0) δ + x v j (cid:1) + M − ∑ j = (cid:26)(cid:20) τ + + ε (cid:16) u n − j (cid:17) (cid:21) v j + ε u n − j v j + ε v j (cid:27) . (2.24)It is easy to check that S ( v ) is strictly convex with the gradient of it denoted as ∇ S ( v ) : =[ ∂ v S ( v ) ,..., ∂ v M S M ( v )] T turning out to be ∂ v j S ( v ) = v j − u nj + u n − j τ − δ x (cid:16) v j + u n − j (cid:17) + (cid:16) v j + u n − j (cid:17) + ε G (cid:16) v j , u n − j (cid:17) . (2.25)By the strict convexity of S ( v ) , we can get the uniqueness of ∇ S ( v ) =
0, which yields theuniqueness of u n + ∈ X M immediately. Thus, the proof is completed. Remark 2.2.
The solvability of the SIFD1 (2.4) can be obtain similarly to the CNFD (2.3)in Lemma 2.3. There exists a unique solution for the SIFD2 due to the fact that it solves alinear system with a strictly diagonally dominant matrix. The solvability and uniquenessfor (2.6) are straightforward since it is explicit.
In this section, we will establish error bounds of the FDTD methods.
Motivated by the analytical results in [16, 18, 25, 26, 31, 35, 38] and references therein, wemake the following assumptions on the exact solution u of the NKGE (2.1) up to the time t = T / ε : ( A ) u ∈ C ([ T / ε ] ; W ∞ p ) ∩ C ([ T / ε ] ; W ∞ ) ∩ C ([ T / ε ] ; W ∞ ) ∩ C ([ T / ε ] ; L ∞ ) , (cid:13)(cid:13)(cid:13)(cid:13) ∂ r + q ∂ t r ∂ x q u ( x , t ) (cid:13)(cid:13)(cid:13)(cid:13) L ∞ .
1, 0 ≤ r ≤
4, 0 ≤ r + q ≤ L ∞ = L ∞ ([ T / ε ] ; L ∞ ) and W m , ∞ p = { u ∈ W m , ∞ | ∂ l ∂ x l u ( a )= ∂ l ∂ x l u ( b ) , 0 ≤ l < m } for m ≥ M = sup ε ∈ ( ] k u ( x , t ) k L ∞ and the grid ‘error’ function e n ∈ X M ( n ≥ ) as e nj = u ( x j , t n ) − u nj , j = M , n = u n ∈ X M is the numerical approximation of the NKGE (2.1).For the CNFD (2.3), we can establish the following error estimates (see its detailedproof in Section 3.2): Theorem 3.1.
Under the assumption (A), there exist constants h > and τ > sufficientlysmall and independent of ε , such that, for any < ε ≤ , when < h ≤ h ε β /2 and < τ ≤ τ ε β /2 ,we have the following error estimates for the CNFD (2.3) with (2.8) and (2.9) k e n k l + k δ + x e n k l . h ε − β + τ ε − β , k u n k l ∞ ≤ + M , 0 ≤ n ≤ T ε − β / τ . (3.2)For the LFFD (2.6), the error estimates can be established as follows (see its detailedproof in Section 3.3): Theorem 3.2.
Assume τ ≤ min { h } and under the assumption (A), there exist constantsh > and τ > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 and < τ ≤ τ ε β /2 and under the stability condition (2.13), we have the errorestimates for the LFFD (2.6) with (2.8) and (2.9) as k e n k l + k δ + x e n k l . h ε − β + τ ε − β , k u n k l ∞ ≤ + M , 0 ≤ n ≤ T ε − β / τ . (3.3)Similarly, for the SIFD1 (2.4) and SIFD2 (2.5), we have the following error estimates(their proofs are quite similar and thus they are omitted for brevity): Theorem 3.3.
Assume τ . h and under the assumption (A), there exist constants h > and τ > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 , < τ ≤ τ ε β /2 and under the stability condition (2.11), we have the following error estimates forthe SIFD1 (2.4) with (2.8) and (2.9) k e n k l + k δ + x e n k l . h ε − β + τ ε − β , k u n k l ∞ ≤ + M , 0 ≤ n ≤ T ε − β / τ . (3.4) Theorem 3.4.
Assume τ . h and under the assumption (A), there exist constants h > and τ > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 , < τ ≤ τ ε β /2 and under the stability condition (2.12), we have the following error estimates forthe SIFD2 (2.5) with (2.8) and (2.9) k e n k l + k δ + x e n k l . h ε − β + τ ε − β , k u n k l ∞ ≤ + M , 0 ≤ n ≤ T ε − β / τ . (3.5) Remark 3.1.
In 2D with d = d = < h . ε β /2 p C d ( h ) and 0 < τ . ε β /2 p C d ( h ) where C d ( h ) = | ln h | when d =
2, and C d ( h ) = h when d = O ( ε − β ) with 0 ≤ β ≤
2. In fact, given an accuracy bound δ >
0, the ε -scalability (or meshing strategy) of the FDTDmethods should be taken as h = O ( ε β /2 p δ ) = O ( ε β /2 ) , τ = O ( ε β /2 p δ ) = O ( ε β /2 ) , 0 < ε ≤
1. (3.6)This implies that, in order to get “correct” numerical solution up to the time at O ( ε − ) ,one has to take the meshing strategy: h = O ( ε ) and τ = O ( ε ) ; and resp., in order toget “correct” numerical solution up to the time at O ( ε − ) , one has to take the meshingstrategy: h = O ( ε ) and τ = O ( ε ) . These results are very useful for practical computationson how to select mesh size and time step such that the numerical results are trustable! For the CNFD (2.3), we establish the error estimates in Theorem 3.1. The key of the proofis to deal with the nonlinearity and overcome the main difficulty in uniformly boundingthe numerical solution u n , i.e., k u n k l ∞ .
1. Here, we adapt the cut-off technique which hasbeen widely used in the literature [1, 2, 39], i.e., the nonlinearity is truncated to a globalLipschitz function with compact support.Denote B = ( + M ) , choose a smooth function ρ ( θ ) ∈ C ∞ ( R + ) and define F B ( θ ) = ρ ( θ / B ) θ , θ ∈ R + , ρ ( θ ) =
1, 0 ≤ θ ≤ ∈ [ ] , 1 ≤ θ ≤ θ ≥
2, (3.7)then F B ( θ ) has compact support and is smooth and global Lipschitz, i.e., there exists C B independent of h , τ and ε , such that | F B ( θ ) − F B ( θ ) | ≤ C B | p θ − p θ | , ∀ θ , θ ∈ R + . (3.8)Set ˆ u = u , ˆ u = u and determine ˆ u n + ∈ X M for n ≥ δ t ˆ u nj − δ x [[ ˆ u ]] nj +[[ ˆ u ]] nj + ε (cid:16) F B (( ˆ u n + j ) )+ F B (( ˆ u n − j ) ) (cid:17) [[ ˆ u ]] nj = j = M −
1. (3.9)In fact, ˆ u nj can be viewed as another approximation of u ( x j , t n ) for j = M and n ≥ τ byusing the properties of ρ and standard techniques in Section 2. Define the corresponding‘error’ function ˆ e n ∈ X M asˆ e nj = u ( x j , t n ) − ˆ u nj , j = M , n ≥
0, (3.10)and we can establish the following estimates: Theorem 3.5.
Under the assumption (A), there exist constants h > and τ > sufficientlysmall and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 and < τ ≤ τ ε β /2 ,we have the following error estimates k ˆ e n k l + k δ + x ˆ e n k l . h ε − β + τ ε − β , k ˆ u n k l ∞ ≤ + M , 0 ≤ n ≤ T ε − β / τ . (3.11)We begin with the local truncation error ˆ ξ n ∈ X M of the scheme (3.9) given asˆ ξ j : = δ + t u ( x j ,0 ) − γ ( x j ) − τ (cid:2) δ x φ ( x j ) − φ ( x j ) − ε ( φ ( x j )) (cid:3) , j = M − ξ nj : = δ t u ( x j , t n ) − (cid:2) δ x u ( x j , t n + )+ δ x u ( x j , t n − ) (cid:3) + (cid:2) u ( x j , t n + )+ u ( x j , t n − ) (cid:3) + ε (cid:0) F B ( u ( x j , t n + ) )+ F B ( u ( x j , t n − ) ) (cid:1)(cid:0) u ( x j , t n + )+ u ( x j , t n − ) (cid:1) , n ≥
1. (3.12)The following estimates hold for ˆ ξ n . Lemma 3.1.
Under the assumption (A), we have k ˆ ξ k l + k δ + x ˆ ξ k l . h + τ , k ˆ ξ n k l . h + τ , 1 ≤ n ≤ T ε − β / τ −
1. (3.13)
Proof.
Under the assumption (A), by applying the Taylor expansion to (3.12), it leads to | ˆ ξ j | . τ k ∂ ttt u k L ∞ + h τ k φ ′′′ k L ∞ . h + τ , j = M − | ˆ ξ nj | . τ (cid:2) k ∂ tttt u k L ∞ + k ∂ ttxx u k L ∞ +( + ε k u k L ∞ ) k ∂ tt u k L ∞ + ε k u k L ∞ k ∂ t u k L ∞ (cid:3) + h k ∂ xxxx u k L ∞ . h + τ , n ≥ | δ + x ˆ ξ j | . h + τ for 0 ≤ j ≤ M −
1. These immediately imply (3.13).Next, we control the nonlinear term as follows.
Lemma 3.2.
For j = M and ≤ n ≤ T ε − β / τ − , denote the error of the nonlinear term ˆ η nj = ε (cid:0) F B ( u ( x j , t n + ) )+ F B ( u ( x j , t n − ) ) (cid:1)(cid:0) u ( x j , t n + )+ u ( x j , t n − ) (cid:1) − ε (cid:16) F B (( ˆ u n + j ) )+ F B (( ˆ u n − j ) ) (cid:17)(cid:16) ˆ u n + j + ˆ u n − j (cid:17) , (3.14) under the assumption (A), we have k ˆ η n k l . ε (cid:16) k ˆ e n − k l + k ˆ e n + k l (cid:17) . (3.15) Proof.
Noticing (3.8) and (3.14), direct calculation for j = M and 1 ≤ n ≤ T ε − β / τ − | ˆ η nj | ≤ C ε h M + | F B (( ˆ u n + j ) ) | + | F B (( ˆ u n − j ) ) | i(cid:16) | ˆ e n + j | + | ˆ e n − j | (cid:17) , (3.16)where the constant C is independent of h , τ and ε . Under the assumption (A) and theproperties of F B , we have k ˆ η n k l . ε h k ˆ e n + k l + k ˆ e n − k l i , 1 ≤ n ≤ T ε − β / τ −
1, (3.17)which completes the proof.Now, we proceed to study the growth of the errors and verify Theorem 3.5. Subtract-ing (3.9) from (3.12), the error ˆ e n ∈ X M satisfies δ t ˆ e nj − (cid:16) δ x ˆ e n + j + δ x ˆ e n − j (cid:17) + (cid:16) ˆ e n + j + ˆ e n − j (cid:17) = ˆ ξ nj − ˆ η nj , 1 ≤ n ≤ T ε − β / τ − e j =
0, ˆ e j = τ ˆ ξ j , j = M −
1. (3.18)Define the ‘energy’ for the error vector ˆ e n asˆ S n = k δ + t ˆ e n k l + (cid:16) k δ + x ˆ e n k l + k δ + x ˆ e n + k l (cid:17) + (cid:16) k ˆ e n k l + k ˆ e n + k l (cid:17) , n ≥
0. (3.19)It is easy to see that ˆ S = k ˆ ξ k l + τ k δ + x ˆ ξ k l + τ k ˆ ξ k l . (cid:0) h + τ (cid:1) . (3.20) Proof. ( Proof of Theorem 3.5 ) When n =
0, the estimates in (3.11) are obvious and the n = < τ < τ and 0 < h < h . Thus,we only need to prove (3.11) for 2 ≤ n ≤ T ε − β / τ .Multiplying both sides of (3.18) by h (cid:16) ˆ e n + j − ˆ e n − j (cid:17) , summing up for j , noticing the fact0 ≤ β ≤ S n − ˆ S n − = h M − ∑ j = (cid:16) ˆ ξ nj − ˆ η nj (cid:17)(cid:16) ˆ e n + j − ˆ e n − j (cid:17) ≤ τε − β (cid:0) k ˆ ξ n k l + k ˆ η n k l (cid:1) + τε β (cid:16) k δ + t ˆ e n k l + k δ + t ˆ e n − k l (cid:17) . ε β τ (cid:16) ˆ S n + ˆ S n − (cid:17) + τε − β (cid:0) h + τ (cid:1) , 1 ≤ n ≤ T ε − β / τ −
1. (3.21)Summing the above inequalities for time steps from 1 to n , there exists a constant C > S n ≤ ˆ S + C ε β τ n ∑ m = ˆ S m + CT ε − β (cid:0) h + τ (cid:1) , 1 ≤ n ≤ T ε − β / τ −
1. (3.22) Hence, the discrete Gronwall’s inequality suggests that there exists a constant τ > < τ ≤ τ , the following holdsˆ S n ≤ (cid:16) ˆ S + CT ε − β (cid:0) h + τ (cid:1) (cid:17) e C ( n + ) ε β τ . ε − β (cid:0) h + τ (cid:1) , 1 ≤ n ≤ T ε − β / τ −
1. (3.23)Recalling k ˆ e n + k l + k δ + x ˆ e n + k l ≤ S n when 0 < ε ≤
1, we can obtain the error estimate k ˆ e n + k l + k δ + x ˆ e n + k l . h ε − β + τ ε − β , 1 ≤ n ≤ T ε − β / τ −
1. (3.24)Finally, we estimate k ˆ u n + k l ∞ for 1 ≤ n ≤ T ε − β / τ −
1. The discrete Sobolev inequalityimplies k ˆ e n k l ∞ ≤ k ˆ e n k l + k δ + x ˆ e n k l . h ε − β + τ ε − β . (3.25)Thus, there exist h > τ > < h ≤ h ε β /2 and 0 < τ ≤ τ ε β /2 ,we obtain k ˆ u n k l ∞ ≤ k u ( x , t n ) k L ∞ + k ˆ e n k l ∞ ≤ M +
1. (3.26)The proof is completed by choosing h = min { h , h } and τ = min { τ , τ , τ } . Proof. ( Proof of Theorem 3.1 ) In view of the definition of ρ , Theorem 3.5 implies that(3.9) collapses to (2.3). By the unique solvability of the CNFD, ˆ u n is identical to u n . Thus,Theorem 3.1 is a direct consequence of Theorem 3.5. For the LFFD (2.6), we establish the error estimates in Theorem 3.2. Throughout thissection, the stability condition (2.13) is assumed. Here, we sketch the proof and omitthose parts similar to the proof of Theorem 3.1 in Section 3.2.
Proof.
Denote the local truncation error as ˜ ξ n ∈ X M ˜ ξ j : = δ + t u ( x j ,0 ) − γ ( x j ) − τ (cid:2) δ x φ ( x j ) − φ ( x j ) − ε φ ( x j ) (cid:3) , j = M − ξ nj : = δ t u ( x j , t n ) − δ x u ( x j , t n )+ u ( x j , t n )+ ε u ( x j , t n ) , 1 ≤ n ≤ T ε − β / τ −
1, (3.27)and the error of the nonlinear term as ˜ η n ∈ X M ˜ η nj : = ε (cid:16) u ( x j , t n ) − ( u nj ) (cid:17) , j = M −
1, 1 ≤ n ≤ T ε − β / τ −
1. (3.28)Similar to Lemma 3.1, under the assumption (A), we have k ˜ ξ k l + k δ + x ˜ ξ k l . h + τ , k ˜ ξ n k l . h + τ , 1 ≤ n ≤ T ε − β / τ −
1. (3.29)The error equation for the LFFD (2.6) can be derived as δ t e nj − δ x e nj + e nj = ˜ ξ nj − ˜ η nj , 1 ≤ n ≤ T ε − β / τ − e j = e j = τ ˜ ξ j , j = M −
1. (3.30) We adapt the mathematical induction to prove Theorem 3.2, i.e. we want to demon-strate that there exist h > τ >
0, such that, when 0 < h < h and 0 < τ < τ , under thestability condition (2.13), the error bounds hold k e n k l + k δ + x e n k l ≤ C (cid:16) h ε − β + τ ε − β (cid:17) , k u n k l ∞ ≤ + M , (3.31)for all 0 ≤ n ≤ T ε − β / τ and 0 ≤ β ≤
2, where C , τ and h will be classified later. For n = n =
1, the error equation (3.30) and the estimate (3.29) imply k e k l = τ k ˜ ξ k l ≤ C τ ( h + τ ) , k δ + x e k l = τ k δ + x ˜ ξ k l ≤ C τ ( h + τ ) . (3.32)In view of the triangle inequality, discrete Sobolev inequality and the assumption (A),there exist h > τ > < h ≤ h and 0 < τ ≤ τ , we have k u k l ∞ ≤ k u ( x , t ) k L ∞ + k e k l ∞ ≤ k u ( x , t ) k L ∞ + k e k l + k δ + x e k l ≤ M +
1. (3.33)In other words, (3.31) hold for n = ≤ n ≤ m − ≤ T ε − β / τ −
1, then we need toshow shat it is still valid when n = m . From (3.28), the error of the nonlinear term can becontrolled as k ˜ η n k l ≤ C ε k e n k l , 1 ≤ n ≤ m −
1. (3.34)Define the ‘energy’ for the error vector e n ( n = ) as S n : = (cid:18) − τ − τ h (cid:19) k δ + t e n k l + n + ∑ k = n k e k k l + h M − ∑ j = (cid:20)(cid:16) e n + j + − e nj (cid:17) + (cid:16) e nj + − e n + j (cid:17) (cid:21) ,where S = (cid:18) − τ − τ h (cid:19) k δ + t e k l + (cid:18) + h (cid:19) k e k l = k ˜ ξ k l ≤ C ( τ + h ) .Under the assumption τ ≤ min { h } , we have 1 − τ /2 − τ / h ≥ >
0. Since k δ + x e n + k l = h M − ∑ j = ( e n + j + − e nj − τδ + t e nj ) ≤ h M − ∑ j = ( e n + j + − e nj ) + τ h k δ + t e n k l ,we can conclude that S n ≥ k δ + x e n + k l + (cid:16) k e n k l + k e n + k l (cid:17) , 1 ≤ n ≤ m −
1. (3.35)Similar to the proof in Section 3.2, there exists τ > < τ ≤ τ , S n ≤ C (cid:16) h ε − β + τ ε − β (cid:17) , 1 ≤ n ≤ m −
1, (3.36) where C depends on T and the exact solution u ( x , t ) . Letting n = m , we have k e m k l + k δ + x e m k l ≤ C ( h ε − β + τ ε − β ) , 1 ≤ m ≤ T ε − β / τ (3.37)where C depends on T and the exact solution u ( x , t ) .It remains to estimate k u m k l ∞ for n = m . In fact, the discrete Sobolev inequality implies k e m k l ∞ ≤ k e m k l + k δ + x e m k l . h ε − β + τ ε − β . (3.38)Thus, there exist h > τ > < h ≤ h ε β /2 and 0 < τ ≤ τ ε β /2 ,we obtain k u m k l ∞ ≤ k u ( x , t m ) k L ∞ + k e m k l ∞ ≤ M +
1, 1 ≤ m ≤ T ε − β / τ . (3.39)Under the stability condition (2.13) and the choices of h = min { h , h } , τ = min { τ , τ , τ } and C = max { C , C } , the estimates in (3.31) are valid when n = m . Hence, the mathemat-ical induction process is done and the proof of Theorem 3.2 is completed. In this section, we present numerical results of the FDTD methods for the NKGE (2.1)up to the long time at O ( ε − β ) with 0 ≤ β ≤ a = b = π and choose the initial data as φ ( x ) = cos ( x )+ cos ( x ) , γ ( x ) = sin ( x ) , 0 ≤ x ≤ π . (4.1)The ‘exact’ solution is obtained numerically by the exponential-wave integrator Fourierpseudospectral method [5, 19] with a very fine mesh size and a very small time step, e.g. h e = π /2 and τ e = − . Denote u nh , τ as the numerical solution at time t = t n obtained bya numerical method with mesh size h and time step τ . In order to quantify the numericalresults, we define the error function as follows: e h , τ ( t n ) = q k u ( · , t n ) − u nh , τ k l + k δ + x ( u ( · , t n ) − u nh , τ ) k l . (4.2)Here we study the following three cases with respect to different 0 ≤ β ≤ O ( ) , i.e., β = O ( ε − ) , i.e., β = O ( ε − ) , i.e., β = t ε = ε β for different 0 < ε ≤
1. In orderto do this, we fix the time step as τ e = − such that the temporal error can be ignored,and solve the NKGE (2.1) under different mesh size h . Tables 1, 3 and 5 depict the spatialerrors for β = β = β =
2, respectively. Then we check the temporal errors at t ε = ε β Table 1: Spatial errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h , τ e ( t = ) h = π /16 h /2 h /2 h /2 h /2 h /2 ε = ε /2 3.33E-2 8.35E-3 2.09E-3 5.22E-4 1.31E-4 3.34E-5order - 2.00 2.00 2.00 1.99 1.97 ε /2 ε /2 ε /2 Table 2: Temporal errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h e , τ ( t = ) τ = τ /2 τ /2 τ /2 τ /2 τ /2 ε = ε /2 2.10E-2 5.45E-3 1.39E-3 3.49E-4 8.76E-5 2.20E-5order - 1.96 1.97 1.99 1.99 1.99 ε /2 ε /2 ε /2 Table 3: Spatial errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h , τ e ( t = ε ) h = π /16 h /2 h /2 h /2 h /2 h /2 ε = - ε /4 7.31E-2 ε /4 ε /4 ε /4 Table 4: Temporal errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h e , τ ( t = ε ) τ = τ /2 τ /2 τ /2 τ /2 τ /2 ε = - ε /4 4.01E-2 ε /4 ε /4 ε /4 Table 5: Spatial errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h , τ e ( t = ε ) h = π /16 h /2 h /2 h /2 h /2 h /2 ε = - ε /2 3.98E-2 ε /2 ε /2 ε /2 Table 6: Temporal errors of the CNFD (2.3) for the NKGE (2.1) with a = , b = π , β = and (4.1) e h e , τ ( t = ε ) τ = τ /2 τ /2 τ /2 τ /2 τ /2 ε = - ε /2 2.56E-2 ε /2 ε /2 ε /2 for different 0 < ε ≤ τ and a fine mesh size h e = π /2 such thatthe spatial errors can be neglected. Tables 2, 4 and 6 show the temporal errors for β = β = β =
2, respectively.From Tables 1-6 for the CNFD and additional similar numerical results for otherFDTD methods not shown here for brevity, we can draw the following observations:(i) For any fixed ε = ε > β =
0, the FDTD methods are uniformly second-orderaccurate in both spatial and temporal discretizations (cf. Tables 1 & 2 and the first rowsin Tables 3-6), which agree with those results in the literature. (ii) In the intermediatelong time regime, i.e. β =
1, the second order convergence in space and time of the FDTDmethods can be observed only when 0 < h . ε and 0 < τ . ε (cf. upper triangles abovethe diagonals (corresponding to h ∼ ε and τ ∼ ε , and being labelled in bold letters)in Tables 3-4), which confirm our error bounds. (iii) In the long time regime, i.e. β =
2, thesecond order convergence in space and time of the FDTD methods can be observed onlywhen 0 < h . ε and 0 < τ . ε (cf. upper triangles above the diagonals (corresponding to h ∼ ε and τ ∼ ε , and being labelled in bold letters) in Tables 5-6), which again confirm ourerror bounds. In summary, our numerical results confirm our rigorous error bounds andshow that they are sharp. Introducing a rescaling in time by s = ε β t with 0 ≤ β ≤ v ( x , s ) : = u ( x , s / ε β )= u ( x , t ) , we can reformulate the NKGE (1.1) into the following oscillatory NKGE ε β ∂ ss v ( x , s ) − ∆ v ( x , s )+ v ( x , s )+ ε v ( x , s ) = x ∈ T d , s > v ( x ,0 ) = φ ( x ) , ∂ s v ( x ,0 ) = ε − β γ ( x ) , x ∈ T d . (5.1)Again, the oscillatory NKGE (5.1) is time symmetric or time reversible and conserves theenergy [5, 19], i.e., E ( s ) : = Z T d (cid:20) ε β | ∂ s v ( x , s ) | + |∇ v ( x , s ) | + | v ( x , s ) | + ε | v ( x , s ) | (cid:21) d x ≡ Z T d (cid:20) | γ ( x ) | + |∇ φ ( x ) | + | φ ( x ) | + ε | φ ( x ) | (cid:21) d x = E ( ) = O ( ) , s ≥
0. (5.2)In fact, the long time dynamics of the NKGE (1.1) up to the time at t = O ( ε − β ) is equivalentto the dynamics of the oscillatory NKGE (5.1) up to the fixed time at s = O ( ) . Of course,the solution of of the NKGE (1.1) propagates waves with wavelength at O ( ) in bothspace and time, and wave speed in space at O ( ) too. On the contrary, the solutionof the oscillatory NKGE (5.1) propagates waves with wavelength at O ( ) in space and O ( ε β ) in time, and wave speed in space at O ( ε − β ) . To illustrate this, Figures 1 & 2 showthe solutions v ( s ) and v ( x ,1 ) , respectively, of the oscillatory NKGE (5.1) with d = T = ( π ) and initial data (4.1) for different 0 < ε ≤ β . We remark here that the Figure 1: The solution v ( s ) of the oscillatory NKGE (5.1) with d = and initial data (4.1) for different ε and β : (a) β = , (b) β = . Figure 2: The solution v ( x ,1 ) of the oscillatory NKGE (5.1) with d = and initial data (4.1) for different ε and β : (a) β = , (b) β = . oscillatory nature of the oscillatory NKGE (5.1) is quite different with that of the NKGEin the nonrelativistic limit regime. In fact, in the nonrelativistic limit regime of the NKGE[3–5, 7], the solution propagates waves with wavelength at O ( ) in space and O ( ε ) intime, and wave speed in space at O ( ) !In the following, we extend the FDTD methods and their error bounds for the NKGE(1.1) in previous sections to the oscillatory NKGE (5.1). Again, for simplicity of notations,the FDTD methods and their error bounds are only presented in 1D, and the results can beeasily generalized to high dimensions with minor modifications. In addition, the proofsfor the error bounds are quite similar to those in Sections 2&3, and thus they are omit-ted for brevity. We adopt similar notations as those used in Sections 2&3 except statedotherwise. In 1D, consider the following oscillatory NKGE ε β ∂ ss v ( x , s ) − ∂ xx v ( x , s )+ v ( x , s )+ ε v ( x , s ) = x ∈ Ω = ( a , b ) , s > v ( x ,0 ) = φ ( x ) , ∂ s v ( x ,0 ) = ε − β γ ( x ) , x ∈ Ω = [ a , b ] , (5.3)with periodic boundary conditions. Choose the temporal step size k : = ∆ s > s n : = nk for n ≥
0. Let v nj be the numerical approximation of v ( x j , s n ) for j = M and n ≥
0, and denote thenumerical solution at time s = s n as v n . Introduce the temporal finite difference operatorsas δ + s v nj = v n + j − v nj k , δ − s v nj = v nj − v n − j k , δ s v nj = v n + j − v nj + v n − j k .We consider the following four FDTD methods:I. The Crank-Nicolson finite difference (CNFD) method ε β δ s v nj − δ x (cid:16) v n + j + v n − j (cid:17) + (cid:16) v n + j + v n − j (cid:17) + ε G (cid:16) v n + j , v n − j (cid:17) =
0, 0 ≤ j ≤ M −
1; (5.4)II. A semi-implicit energy conservative finite difference (SIFD1) method ε β δ s v nj − δ x v nj + (cid:16) v n + j + v n − j (cid:17) + ε G (cid:16) v n + j , v n − j (cid:17) =
0, 0 ≤ j ≤ M −
1; (5.5)III. Another semi-implicit finite difference (SIFD2) method ε β δ s v nj − δ x (cid:16) v n + j + v n − j (cid:17) + (cid:16) v n + j + v n − j (cid:17) + ε ( v nj ) =
0, 0 ≤ j ≤ M −
1; (5.6)IV. The Leap-frog finite difference (LFFD) method ε β δ s v nj − δ x v nj + v nj + ε ( v nj ) =
0, 0 ≤ j ≤ M − n ≥
1. (5.7)The initial and boundary conditions are discretized as v n + = v n + M , v n + − = v n + M − , n ≥ v j = φ ( x j ) , j = M . (5.8)Using the Taylor expansion and noticing (5.3), the first step v ∈ X M can be computed as v j = φ ( x j )+ k ε − β γ ( x j )+ k ε − β (cid:2) δ x φ ( x j ) − φ ( x j ) − ε φ ( x j ) (cid:3) , 0 ≤ j ≤ M −
1. (5.9)In fact, if we take k = τε β in the FDTD methods in this section, then they are consistentwith those FDTD methods presented in Section 2. Thus they have the same solutions.We remark here that, in practical computations, in order to uniformly bound the firststep value v ∈ X M for ε ∈ ( ] , in the above approximation (5.9), k ε − β and k ε − β arereplaced by sin ( k ε − β ) and k sin ( k ε − β ) , respectively [5, 8]. Denote ˜ σ max : = max ≤ n ≤ T / k k v n k l ∞ . (5.10)Similar to Section 2, following the von Neumann linear stability analysis of the classicalFDTD methods for the NKGE in the nonrelativistic limit regime [5, 29], we can concludethe linear stability of the above FDTD methods for oscillatory NKGE (5.3) up to the fixedtime s = T in the following lemma. Lemma 5.1.
For the above FDTD methods applied to the oscillatory NKGE (5.3) up to the fixedtime s = T , we have:(i) The CNFD (5.4) is unconditionally stable for any h > k > and < ε ≤ .(ii) When h ≥ , the SIFD1 (5.5) is unconditionally stable for any h > and k > ; and when < h < , this scheme is conditionally stable under the stability condition < k < ε β h / p − h , h >
0, 0 < ε ≤
1. (5.11) (iii) When ˜ σ max ≤ ε − , the SIFD2 (5.6) is unconditionally stable for any h > and k > ; andwhen ˜ σ max > ε − , this scheme is conditionally stable under the stability condition < k < ε β / p ε ˜ σ max − h >
0, 0 < ε ≤
1. (5.12) (iv) The LFFD (5.7) is conditionally stable under the stability condition < k < ε β h / q + h ( + ε ˜ σ max ) , h >
0, 0 < ε ≤
1. (5.13)For the CNFD (5.4) and SIFD1 (5.5), we have the following energy conservation prop-erties:
Lemma 5.2.
The CNFD (5.4) conserves the discrete energy as E n = ε β k δ + s v n k l + (cid:16) k δ + x v n k l + k δ + x v n + k l (cid:17) + (cid:16) k v n k l + k v n + k l (cid:17) + h ε M − ∑ j = h | v nj | + | v n + j | i ≡ E , n = Similarly, the SIFD1 (5.5) conserves the discrete energy as ˜ E n = ε β k δ + s v n k l + h M − ∑ j = (cid:16) δ + x v nj (cid:17)(cid:16) δ + x v n + j (cid:17) + (cid:16) k v n k l + k v n + k l (cid:17) + h ε M − ∑ j = h | v nj | + | v n + j | i ≡ ˜ E , n = Again, motivated by the analytical results and the assumptions on the NKGE (2.1), weassume that the exact solution v of the oscillatory NKGE (5.3) satisfies ( B ) v ∈ C ([ T ] ; W ∞ p ) ∩ C ([ T ] ; W ∞ ) ∩ C ([ T ] ; W ∞ ) ∩ C ([ T ] ; L ∞ ) , (cid:13)(cid:13)(cid:13)(cid:13) ∂ r + q ∂ s r ∂ x q v ( x , s ) (cid:13)(cid:13)(cid:13)(cid:13) L ∞ ([ T ] ; L ∞ ) . ε β r , 0 ≤ r ≤
4, 0 ≤ r + q ≤ e n ∈ X M ( n ≥ ) as˜ e nj = v ( x j , s n ) − v nj , j = M , n = v n ∈ X M is the numerical approximation of the oscillatory NKGE (5.3) obtained byone of the FDTD methods.By taking k = τε β in the above FDTD methods and noting the error bounds in Section3, we can immediately obtain error bounds of the above FDTD methods for the oscillatoryNKGE (5.3). Theorem 5.1.
Under the assumption (B), there exist constants h > and k > sufficientlysmall and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 and < k ≤ k ε β /2 ,we have the following error estimates for the CNFD (5.4) with (5.8) and (5.9) k ˜ e n k l + k δ + x ˜ e n k l . h ε − β + k ε − β , k v n k l ∞ ≤ + M , 0 ≤ n ≤ T / k . (5.17) Theorem 5.2.
Assume k . h ε β and under the assumption (B), there exist constants h > andk > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 , < k ≤ k ε β /2 and under the stability condition (5.11), we have the following error estimates forthe SIFD1 (5.5) with (5.8) and (5.9) k ˜ e n k l + k δ + x ˜ e n k l . h ε − β + k ε − β , k v n k l ∞ ≤ + M , 0 ≤ n ≤ T / k . (5.18) Theorem 5.3.
Assume k . h ε β and under the assumption (B), there exist constants h > andk > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 , < k ≤ k ε β /2 and under the stability condition (5.12), we have the following error estimates forthe SIFD2 (5.6) with (5.8) and (5.9) k ˜ e n k l + k δ + x ˜ e n k l . h ε − β + k ε − β , k v n k l ∞ ≤ + M , 0 ≤ n ≤ T / k . (5.19) Theorem 5.4.
Assume k . h ε β and under the assumption (B), there exist constants h > andk > sufficiently small and independent of ε , such that for any < ε ≤ , when < h ≤ h ε β /2 , < k ≤ k ε β /2 and under the stability condition (5.13), we have the following error estimates forthe LFFD (5.7) with (5.8) and (5.9) k ˜ e n k l + k δ + x ˜ e n k l . h ε − β + k ε − β , k v n k l ∞ ≤ + M , 0 ≤ n ≤ T / k . (5.20)The above four FDTD methods share the same spatial/temporal resolution capacityfor the oscillatory NKGE (5.3) up to the fixed time at O ( ) . In fact, given an accuracybound δ >
0, the ε -scalability (or meshing strategy) of the FDTD methods for the oscilla-tory NKGE (5.3) should be taken as h = O ( ε β /2 p δ ) = O ( ε β /2 ) , τ = O ( ε β /2 p δ ) = O ( ε β /2 ) , 0 < ε ≤
1. (5.21)Again, these results are very useful for practical computations on how to select mesh sizeand time step such that the numerical results are trustable!
Consider the following oscillatory NKGE in d -dimensional ( d = ε β ∂ ss v ( x , s ) − ∆ v ( x , s )+ v ( x , s )+ ε v ( x , s ) = x ∈ R d , s > v ( x ,0 ) = φ ( x ) , ∂ s v ( x ,0 ) = ε − β γ ( x ) , x ∈ R d . (5.22)Similar to the oscillatory NKGE (5.1), the solution of of the oscillatory NKGE (5.22) prop-agates waves with wavelength at O ( ) in space and O ( ε β ) in time, and wave speed inspace at O ( ε − β ) . To illustrate the rapid wave propagation in space at O ( ε − β ) , Figure 3shows the solution v ( x ,1 ) of the oscillatory NKGE (5.22) with d = φ ( x ) = ( e x + e − x ) and γ ( x ) = x ∈ R . (5.23)Similar to those in the literature, by using the fast decay of the solution of the os-cillatory NKGE (5.22) at the far field (see [5, 20, 37] and references therein), in practicalcomputation, we usually truncate the originally whole space problem onto a boundeddomain Ω with periodic boundary conditions, provided that Ω is large enough such thatthe truncation error is negligible. Then the truncated problem can be solved by the FDTDmethods. Of course, due to the rapid wave propagation in space of the oscillatory NKGE(5.22) (cf. Fig. 3), in order to compute numerical solution up to the time at O ( ) , ingeneral, the size of the bounded domain Ω has to be taken as O ( ε − β ) .In the following, we report numerical results of the oscillatory NKGE (5.22) with d = -10 -5 0 5 10-0.8-0.6-0.4-0.200.20.4 -60 -40 -20 0 20 40 60-0.5-0.4-0.3-0.2-0.100.10.20.3 Figure 3: The solutions v ( x ,1 ) of the oscillatory NKGE (5.22) with d = and initial data (5.23) for different ε and β : (a) β = , (b) β = . Ω ε = [ − − ε − β ,4 + ε − β ] . The ‘exact’ solution is obtained numerically by the exponential-wave integrator Fourier pseudospectral method with a very fine mesh size and a verysmall time step, e.g. h e = and k e = × − . Denote v nh , k as the numerical solutionat s = s n obtained by a numerical method with mesh size h and time step k . In order toquantify the numerical results, we define the error function as follows: e h , k ( s n ) = q k v ( · , s n ) − v nh , k k l + k δ + x ( v ( · , s n ) − v nh , k ) k l . (5.24)Tables 7 and 8 show the spatial and temporal errors, respectively, of the CNFD methodwith β =
1, and Tables 9 and 10 show similar results for β =
2. The results for other FDTDmethods are quite similar and they are omitted here for brevity.From Tables 7-10 for the CNFD and additional similar numerical results for otherFDTD methods not shown here for brevity, we can draw the following observations onthe FDTD methods for the oscillatory NKGE (5.1) (or (5.22)):(i) For any fixed ε = ε > β =
0, the FDTD methods are uniformly second-orderaccurate in both spatial and temporal discretizations (cf. the first rows in Tables 7-10),which agree with those results in the literature. (ii) In the intermediate oscillatory case, Table 7: Spatial errors of the CNFD (5.4) for the oscillatory NKGE (5.22) with d = , β = and (5.23) e h , k e ( s = ) h = h /2 h /2 h /2 h /2 h /2 ε = - ε /4 5.60E-2 ε /4 ε /4 ε /4 Table 8: Temporal errors of the CNFD (5.4) for the oscillatory NKGE (5.22) with d = , β = and (5.23) e h e , k ( s = ) k = k /4 k /4 k /4 k /4 k /4 ε = < - ε /4 ε /4 ε /4 ε /4 Table 9: Spatial errors of the CNFD (5.4) for the oscillatory NKGE (5.22) with d = , β = and (5.23) e h , k e ( s = ) h = h /2 h /2 h /2 h /2 h /2 ε = - ε /2 5.64E-2 ε /2 ε /2 ε /2 Table 10: Temporal errors of the CNFD (5.4) for the oscillatory NKGE (5.22) with d = , β = and (5.23) e h e , k ( s = ) k = k /4 k /4 k /4 k /4 k /4 ε = < - ε /4 ε /4 ε /4 ε /4 β =
1, the second order convergence in space and time of the FDTD methods can beobserved only when 0 < h . ε and 0 < k . ε (cf. upper triangles above the diagonals(corresponding to h ∼ ε and k ∼ ε , and being labelled in bold letters) in Tables 7-8),which confirm our error bounds. (iii) In the highly oscillatory case, i.e. β =
2, the secondorder convergence in space and time of the FDTD methods can be observed only when0 < h . ε and 0 < k . ε (cf. upper triangles above the diagonals (corresponding to h ∼ ε and k ∼ ε , and being labelled in bold letters) in Tables 9-10), which again confirm our errorbounds. In summary, our numerical results confirm our rigorous error bounds and showthat they are sharp. Four different finite difference time domain FDTD methods were adapted to discretizethe nonlinear Klein-Gordon equation (NKGE) with a weak cubic nonlinearity, while thenonlinearity strength is characterized by ε with 0 < ε ≤ O ( ε − β ) with 0 ≤ β ≤
2. The error bounds depend explicitly on the mesh size h andtime step τ as well as the small parameter ε ∈ ( ] , which indicate the temporal and spa-tial resolution capacities of the FDTD methods for the long time dynamics of the NKGE.Based on the error bounds, in order to get “correct” numerical solution of the NKGE upto the long time at O ( ε − β ) with 0 < β ≤
2, the ε -scalability (or meshing strategy) of theFDTD methods has to be taken as: h = O ( ε β /2 ) and τ = O ( ε β /2 ) . In addition, the FDTDmethods were also applied to solve an oscillatory NKGE and their error bounds werealso obtained. Extensive numerical results were reported to confirm our error boundsand to demonstrate that they are sharp. Acknowledgments
We thank fruitful discussion with Dr Chunmei Su. This work was partially supported bythe Ministry of Education of Singapore grant R-146-000-223-112.
References [1] W. B
AO AND
Y. C AI , Uniform error estimates of finite difference methods for the nonlinearSchr ¨odinger equation with wave operator, SIAM J. Numer. Anal., 50 (2012) 492-521.[2] W. B AO AND
Y. C AI , Optimal error estimates of finite difference methods for the Gross-Pitaevskii equation with angular momentum rotation, Math. Comp., 82 (2013) 99-128.[3] W. B AO , Y. C AI , X. J IA AND
J. Y IN , Error estimates of numerical methods for the nonlinearDirac equation in the nonrelativistic limit regime, Sci. China Math., 59 (2016), 1461-1494.[4] W. B AO , Y. C AI AND
X. Z
HAO , A uniformly accurate multiscale time integrator pseudospec-tral method for the Klein-Gordon equation in the nonrelativistic limit regime, SIAM J. Nu-mer. Anal., 52 (2014) 2488-2511.[5] W. B
AO AND
X. D
ONG , Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime, Numer. Math., 120 (2012) 189-229.[6] W. B AO , X. D ONG AND
X. Z
HAO , An exponential wave integrator pseudospectral methodfor the Klein-Gordon-Zakharov system, SIAM J. Sci. Comput., 35 (2013), A2903-A2927.[7] W. B AO , X. D ONG AND
X. Z
HAO , Uniformly accurate multiscale time integrators for highlyoscillatory second order differential equations, J. Math. Study, 47 (2014) 111-150.[8] W. B
AO AND
C. S U , Uniform error bounds of a finite difference method for the Klein-Gordon-Zakharov system in the subsonic limit regime, Math. Comp., 87 (2018) 2133-2158.[9] W. B AO AND
L. Y
ANG , Efficient and accurate numerical methods for the Klein-Gordon-Schr ¨odinger equations, J. Comput. Phys., 225 (2007), 1863-1893.[10] J. B
OURGAIN , Construction of approximative and almost periodic solutions of perturbedlinear Schr ¨odinger and wave equations, Geom. Funct. Anal., 6 (1996) 201-230.[11] P. B
RENNER , On the existence of global smooth solutions of certain semi-linear hyperbolicequations, Math. Z., 167 (1979) 99-135.[12] P. B
RENNER AND W. VON W AHL , Global classical solutions of nonlinear Klein-Gordon equa-tions, Math. Z., 176 (1981) 87-121.[13] F.E. B
ROWDER , On nonlinear Klein-Gordon equations, Math. Z., 80 (1962) 249-264.[14] Q. C
HANG , G. W
ANG AND
B. G UO , Conservative scheme for a model of nonlinear disper-sive waves and its solitary waves induced by boundary motion, J. Comput. Phys., 93 (1991)360-375.[15] S. C. C HIKWENDU AND
C. V. E
ASWARAN , Multiple-scale solution of initial-boundary valueproblems for weakly nonlinear Klein-Gordon equations on the semi-infinite line, SIAM J.Appl. Math., 52 (1992) 946-958.[16] J.-M. D
ELORT , Temps d’existence pour l’´equation de Klein-Gordon semi-lin´eaire `a donn´eespetites p´eriodiques, Amer. J. Math., 120 (1998) 663-689.[17] J.-M. D
ELORT , On long time existence for small solutions of semi-linear Klein-Gordon equa-tions on the torus, J. Anal. Math., 107 (2009) 161-194.[18] J.-M. D
ELORT AND
J. S
ZEFTEL , Long time existence for small data nonlinear Klein-Gordonequations on tori and spheres, Int. Math. Res. Not. IMRN, 37 (2004) 1897-1966. [19] X. D ONG , Z. X
U AND
X. Z
HAO , On time-splitting pseudospectral discretization for non-linear Klein-Gordon equation in nonrelativistic limit regime, Commun. Comput. Phys., 16(2014) 440-466.[20] D. B. D
UNCAN , Sympletic finite difference approximations of the nonlinear Klein-Gordonequation, SIAM J. Numer. Anal., 34 (1997) 1742-1760.[21] D. F
ANG AND
Q. Z
HANG , Long-time existence for semi-linear Klein–Gordon equations ontori, J. Differential Equations, 249 (2010) 151-179.[22] E. F
AOU AND
K. S
CHRATZ , Asymptotic preserving schemes for the Klein–Gordon equationin the nonrelativistic limit regime, Numer. Math., 126 (2014) 441-469.[23] H. F
ESHBACH AND
F. V
ILLARS , Elementary relativistic wave mechanics of spin 0 and spin1/2 particles, Rev. Modern Phys., 30 (1958) 24.[24] S. J IM ´ ENEZ AND
L. V ´
AZQUEZ , Analysis of four numerical schemes for a nonlinear Klein-Gordon equation, Appl. Math. Comput., 35 (1990) 61-94.[25] M. K
EEL AND
T. T AO , Small data blow-up for semilinear Klein-Gordon equations, Amer. J.Math., 121 (1999) 629-669.[26] S. K LAINERMAN , Global existence of small amplitude solutions to nonlinear Klein-Gordonequations in four space-time dimensions, Comm. Pure Appl. Math., 38 (1985) 631-641.[27] P. S. L
ANDA , Nonlinear oscillations and waves in dynamical systems, Kluwer AcademicPublishers, Boston, MA, 1996.[28] R. L
ANDES , On Galerkin’s method in the existence theory of quasilinear elliptic equations,J. Funct. Anal., 39 (1980) 123-148.[29] R. J. L
EVEQUE , Finite volume methods for hyperbolic problems, Cambridge UniversityPress, the United Kingdom, 2002.[30] S. L
I AND
L. V U -Q UOC , Finite difference calculus invariant structure of a class of algorithmsfor the nonlinear Klein-Gordon equation, SIAM J. Numer. Anal., 32 (1995) 1839-1875.[31] T. L
I AND
Y. Z
HOU , Nonlinear Klein-Gordon equations, Springer, Berlin, 2017.[32] H. L
INDBLAD , On the lifespan of solutions of nonlinear Klein-Gordon equations with smallinitial data, Comm. Pure Appl. Math., 43 (1990) 445-472.[33] S. M
ACHIHARA , The nonrelativistic limit of the nonlinear Klein-Gordon equation, Funkcial.Ekvac., 44 (2001) 243-252.[34] S. M
ACHIHARA , K. N
AKANISHI AND
T. O
ZAWA , Nonrelativistic limit in the energy spacefor nonlinear Klein-Gordon equations, Math. Ann., 322 (2002) 603-621.[35] K. O NO , Global existence and asymptotic behavior of small solutions for semilinear dissi-pative Klein-Gordon equations, Discrete Contin. Dyn. Syst., 9 (2003) 651-662.[36] J. J. S AKURAI , Advanced Quantum Mechanics, Addison-Wesley, New York, 1967.[37] W. S
TRAUSS AND
L. V ´
AZQUEZ , Numerical solution of a nonlinear Klein-Gordon equation,J. Comput. Phys., 28 (1978) 271-278.[38] G. T
ODOROVA AND
B. Y
ORDANOV , Critical exponent for a nonlinear Klein-Gordon equa-tion with damping, J. Differential Equations, 174 (2001) 464-489.[39] V. T
HOM ´ EE , Galerkin finite element methods for parabolic problems, Springer, Berlin, 1997.[40] W. T. V AN H ORSSEN , An asymptotic theory for a class of initial-boundary value problemsfor weakly nonlinear wave equations with an application to a model of the galloping oscil-lations of overhead transmission lines, SIAM J. Appl. Math., 48 (1988) 1227-1243.[41] F. V
ERHULST , Methods and Applications of Singular Perturbations: Boundary Layers andMultiple Timescale Dynamics, Texts Appl. Math. 50, Springer, New York, 2005.[42] W.
VON W AHL , Regular solutions of initial-boundary value problems for linear and nonlin-ear wave-equations. II, Math. Z., 142 (1975) 121-130. [43] D. W ILLETT AND
J. W
ONG , On the discrete analogues of some generalizations of Gronwall’sinequality, Monatsh. Math., 69 (1965) 362-367.[44] L. Z
HANG , Convergence of a conservative difference scheme for a class of Klein–Gordon–Schr ¨odinger equations in one space dimension, Appl. Math. Comput., 163 (2005) 343-355.[45] Y. Z