- Email: [email protected]
Large amplitude instability in finite difference approximations to the Klein–Gordon equation
Applied Numerical Mathematics 31 (1999) 173–182
Large amplitude instability in finite difference approximations to the Klein–Gordon equation Mark A.M. Lynch 1 Glasgow Caledonian University, Mathematics Department, City Campus, Cowcaddens Road, Glasgow G4 0BA, Scotland, UK
Abstract In a paper by Jiminez and Vazquez (1990), four finite difference schemes for approximating the nonlinear Klein– Gordon equation were discussed. They observed undesirable characteristics in some of the numerical schemes, in particular a loss of spatial symmetry and the onset of instability for large values of a parameter in the initial condition of the equation. The proposition, put forward by Jiminez and Vazquez, that the observed instability is protected against by using energy conserving schemes is questioned. An analysis of the schemes as applied to a linear problem is carried out and these indicate that the instability arises from the use of explicit finite difference schemes rather than any failure of energy conservation. This conjecture is further supported by an analysis of two further schemes. Some general results applicable to the nonlinear finite difference schemes are also presented. 1999 Elsevier Science B.V. and IMACS. All rights reserved. Keywords: Finite difference; Klein–Gordon; Stability
1. Introduction The numerical solutions to the NonLinear Klein–Gordon Equation (NLKGE) Ut t − Uxx + f (U ) = 0
(1.1)
have received considerable attention in the literature. A variety of second order finite difference schemes have been presented [8,12,13], and a strong emphasis has been given to developing schemes which conserve a discrete energy [7,8]. The latter approach has the particular appeal in that it reflects the energy conservation property of the differential equation. We note that alternative approaches using spectral and pseudospectral methods have recently been presented [4,9]. Investigations have been carried out on the particular case of the NLKGE, called the sine-Gordon equation [3,6], and, in particular, a finite element method is described in [2]. 1 E-mail: [email protected].
0168-9274/99/$20.00 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 8 ) 0 0 1 2 8 - 7
174
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
In a paper by Jiminez and Vazquez [8], four finite difference schemes for approximating the nonlinear Klein–Gordon equation were discussed. They observed undesirable characteristics in some of the numerical schemes, in particular a loss of spatial symmetry and the onset of instability for large values of a parameter in the equation’s initial condition. This instability will be referred to as large amplitude instability in this paper. The proposition, put forward by Jiminez and Vazquez [8], that the onset of this type of instability is protected against by using energy conserving schemes is questioned. In this paper the schemes are applied to the Linear Klein–Gordon Equation (LKGE) and a standard von Neumann analysis is performed to gain insight into the source of this large amplitude instability. The analysis suggests the counter proposition that the protection arises through the implicit nature of the schemes. To test this proposition two further schemes are considered: an explicit scheme with a known conserved energy and an implicit scheme with no known conserved energy. The explicit scheme used is an extension of the one described by Fei and Vazquez [6] for investigating the sine-Gordon equation. The von Neumann analysis suggests that the explicit scheme is not protected against the onset of the large amplitude instability, even though it has a conserved quantity, whereas the implicit scheme is protected from this type of instability. This conclusion is supported by numerical experiments with the two schemes. Before proceeding to an investigation of the particular schemes discussed by Jiminez and Vazquez, some general results are presented. These results cover the accuracy of some of the energy conserving schemes discussed in the literature and sufficiency conditions for the existence of an energy type conservation property. 1.1. Background We introduce some basic notation. Define τ as the time increment and χ is the spacing between points on the mesh (assumed uniform). Let U (t, x) be an exact solution to the NLKGE (1.1). Define a discrete function unj which is an approximation to the exact solution of the NLKGE, i.e., unj ≈ U (nτ, j χ). The schemes discussed in this paper will make use of the following central difference operators: δt2 unj = un+1 − 2unj + un−1 and δx2 unj = unj+1 − 2unj + unj−1 . Finally, we introduce the function F (U ), j j a primitive of the function f (U ). In [8] four second order finite difference (FD) schemes were subjected to a variety of numerical experiments. We concentrate on one of these experiments where the large amplitude of the initial condition resulted in three of the schemes failing to successfully propagate the solution whereas an energy conserving scheme succeeded. It is worth noting that none of the schemes exhibited any instability when the initial condition was of a ‘reasonable’ size. The schemes were set to solve the following initial value problem with periodic boundary condition: Ut t − Uxx + U + U 3 = 0 subject to the initial conditions
2π U (0, x) = A 1 + cos x 1.28
(1.2)
,
Ut (0, x) = 0,
defined on the space interval [0, 1.28]. The time step size was 0.01 and the mesh spacing 0.02. The parameter A was the subject of experimental variation. The four schemes used in the experimentation were:
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
175
(A) due to Strauss and Vazquez [10]: n−1 1 2 n 1 2 n F (un+1 j ) − F (uj ) δ u − δ u + = 0, τ2 t j χ2 x j un+1 − un−1 j j
(B) in common use: 1 2 n 1 δt uj − 2 δx2 unj + f unj = 0, 2 τ χ
(C) due to Ablowitz et al. [1]:
unj+1 + unj−1 1 2 n 1 2 n δ u − δ u +f τ2 t j χ2 x j 2
= 0,
(D) due to Jiminez [7]: 1 2 n 1 2 n F (unj+1 ) − F (unj−1 ) δ u − δ u + = 0. τ2 t j χ2 x j unj+1 − unj−1 The conclusions of the paper were that scheme (A) was superior to the other three in that it better preserved the spatial symmetry of the numerical solution and was not prone to the onset of instability for large values of the amplitude parameter A. Unfortunately, when looking for numerical schemes with higher orders of accuracy, it is not clear how to take advantage of this claim unless a conserved discrete energy can be found. In this paper we argue that this line of attack is not of itself a useful criterion for high order schemes. We use an analysis of the linear problem to propose that it is not the existence of the energy conservation property of the scheme that gives it protection from large amplitude instability, rather it is the implicit nature of the scheme.
2. Accuracy and a conserved quantity The energy conserving scheme discussed in [8] belongs to a family of schemes with the following structure: n+1/2
n−1/2
) − N(uj N(uj 1 2 n 1 δt uj − 2 δt2 unj + n+1 2 τ χ uj − un−1 j
)
= 0.
(2.1)
The symbol N represents some combination of the functions f and F with space and time differencing operators. We give two useful results relating to this generic scheme. Let u and v be discrete square summable functions defined for all integers (Z). Then hu, vi =
X
ui vi χ.
i∈Z
We define the following additional time differencing operators: 1t unj = un+1 − un−1 j j ,
n+1/2
Γt /2un = uj
and the following space differencing operators:
n−1/2
− uj
,
δt2 unj = un+1 − 2unj + un−1 j j ,
176
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
1Fx uj = unj+1 − unj , δx2 unj = unj+1 − 2unj + unj−1 , unj+1/2 + unj−1/2 µx uj = , 1x uj = unj+1 − unj−1 . 2 The application of space differencing operators to the functions u and v will be assumed to generate functions that are always square summable. Theorem 1. The numerical schemes (2.1) are second order accurate approximations to the nonlinear Klein–Gordon equation (1.1) if they satisfy n+1/2
n−1/2
) − N(uj ) ∂ (2.2) = 2 F unj + O τ 2 , τ ∂t where (∂/∂t)F (unj ) is interpreted as the time derivative of F evaluated at an exact solution of the partial differential equation at time j τ and at position j χ . N(uj
Proof. The linear part of scheme (2.1) is already a second order accurate approximation to the linear wave equation. Assume the scheme satisfies (2.2) this implies that n+1/2
n+1/2
) − N(uj ) ∂ n = 2f unj uj + O τ 2 . τ ∂t Using the second order central difference approximation to the time derivative of u and multiplying through by τ implies N(uj
n+1/2
N uj
n−1/2
− N uj
= un+1 − un−1 f unj + O τ 3 . j j
Dividing through by the order τ term un+1 − un−1 completes the proof. j j
2
Theorem 2. Solutions to the finite difference schemes (2.1) satisfy the following conservation condition:
n+1/2
1t /2 unj , 1t /2unj + ρ 2 1Fx uj
−1/2
, 1Fx uj
+ τ 2χ
X
N unj = K,
j
where K is a constant and the mesh ratio ρ is defined by ρ = τ/χ . Proof. Premultiply (2.1) by τ 2 χ1t unj and summing over all integer values j gives
n+1/2
1t unj , δt2 unj + 1t unj , ρ 2 δx2 unj + τ 2 1t unj ,
N(uj
n−1/2
) − N(uj
un+1 − un−1 j j
Application of the Lemmas 1–3 in Appendix reduces this to
1t /2
1t /2 unj , 1t /2 unj
+ρ
2
1Fx unj , 1Fx unj
+τ χ 2
X
N
unj
)
= 0.
= 0.
j
Hence the quantity in the square bracket is constant with respect to time. It is noted that for the particular NLKGE (1.2) the term N(u) will be positive definite and the conserved quantity is a true energy whenever ρ < 1 (see also the remarks by Jiminez [7]).
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182 n+1/2
177
n−1/2
Two particular choices of note are: N(unj ) = (F (uj ) + F (uj )) which returns scheme (A) n+1/2 n−1/2 n above, and N(uj ) = 2F ((uj + uj )/2) which returns a scheme described in [12] which is also a generalization of a scheme described in [6] for the sine-Gordon equation. 2 In Section 4 we will make use of the finite difference scheme described in Theorem 3. Theorem 3. Solutions to the finite difference schemes n+1 N(un+2 1 1 2 n+2 j ) − N(uj ) n+3 n+2 n+1 n n+1 u − uj + uj + uj − 2 δx uj + uj + =0 2τ 2 j 2χ un+2 − un+1 j j
satisfy the following conservation condition:
n+3/2
1t /2 uj
n+1/2
, 1t /2uj
n+3/2
+ ρ 2 1Fx uj
n+3/2
, 1Fx uj
+ 2τ 2 χ
X
n+3/2
N uj
=K
j
for some constant K. Proof. Multiply the difference equation through by 2τ 2 χ(un+2 − un+1 j j ), substitute ρ = τ/χ and sum over all integer values of j . Treating each of the three parts under the summation separately. The time differencing operator:
n+3 n+2 un+2 − un+1 − un+2 + un+1 + unj = un+2 − un+2 − un+1 j j , uj j j j j , 1t uj j
n+3/2
n+1/2
= 1t /2 1t /2 uj , 1t /2uj (applying Lemma 2 in Appendix). The space differencing operator:
2 n+2 un+2 − un+1 + δx2 un+1 j j , δx uj j
2 n+2 2 n+1 2 n+2 2 n+1 = un+2 − un+1 − un+1 + un+2 j , δx uj j , δx uj j , δx uj j , δx uj 2 n+2 2 n+1 = un+2 − un+1 j , δx uj j , δx uj
F n+2 F n+1 = − 1Fx un+2 + 1Fx un+1 j , 1x uj j , 1x uj
=
n+3/2 n+3/2 −1t /2 1Fx uj , 1Fx uj .
The nonlinear term is straightforward.
(applying Lemma 3 in Appendix)
2
It is an immediate consequence of Theorems 1 and 2 that the schemes (A)–(D) are second order and that scheme (A) has a conserved energy. The schemes (B)–(D) do not satisfy the conditions of Theorem 2 but this of itself does not mean that they fail to have a conserved quantity. As far as is known these schemes do not have a conserved energy, though scheme (D) (as does (A)) has a conserved momentum [8]. 3. Linear schemes The nonlinear schemes (A)–(D) are applied to the simpler linear Klein–Gordon equation Ut t − Uxx + b2 U = 0,
(3.1)
178
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
where b is a constant. The resulting schemes are then subjected to a von Neumann stability analysis. The resulting schemes are: (A1) (B1)
1 2 n 1 2 n b2 n+1 n−1 δ u − δ u + u + u = 0, t j x j j j τ2 χ2 2 1 2 n 1 δt uj − 2 δx2 unj + b2 unj = 0, 2 τ χ
(C1/D1) (note that (C) and (D) give rise to the same linear scheme): 1 2 n 1 2 n b2 n δ u − δx uj + uj +1 + unj−1 = 0. t j 2 2 τ χ 2
Applying a standard von Neumann stability analysis to the schemes (A1), (B1) and (C1/D1) results in characteristic equations with the following structure: ξ 2 − λξ + 1 = 0, where λ is determined by the scheme. The characteristic equations have two roots which we signify by ξ1 and ξ2 or simply ξ1,2 . Definition. The linear numerical schemes (A1), (B1), (C1/D1) are stable if and only if |ξ1,2 | 6 1 and ξ1 6= ξ2 whenever |ξ1,2 | = 1. The von Neumann analysis returns the following stability conditions on the mesh ratio ρ = τ/χ : • scheme (A1): ρ 2 < 1 + 14 b2 τ 2 , • scheme (B1) and (C1/D1): ρ 2 < 1 − 14 b2 τ 2 . The above results can be used to gain insight into the instability observed by Jiminez and Vazquez [8]. They observed in their numerical experiments that when the parameter A in the initial condition was allowed to become too large then the schemes (B)–(D) became unstable. Whereas scheme (A) remained stable for all observed values of the parameter A. An inspection of the above results shows that the stability of the schemes (B1) and (C1/D1) are adversely dependent on the parameter b2 , in the sense that increasing the value of b2 (with fixed τ ) eventually results in schemes that are unstable. On the other hand, the stability of the scheme (A1) is seen to be independent of the value of b2 . We note in passing that for both the schemes (B1) and (C1/D1) the onset of large amplitude instability can be compensated for by decreasing the time step size τ . The above arguments can be used to derive a possible explanation for the observations made by Jiminez and Vazquez in their numerical experiments. For large values of the parameter A in the initial condition the function f (U ) = U + U 3 will take values in the range 0 to order A3 . Any reasonable linear approximation of the form f (U ) ≈ b2 U will be a very steeply increasing function. The results of the preceding analysis indicate that when the schemes (B1) and (C1/D1) are applied to a linear function with sufficiently steep gradient numerical instability can be expected. On the other hand, no matter how large the value of b2 , the scheme (A1) will remain stable. This is not to say that other, nonlinear sources of instability, may be forced by the FD schemes but we have gained a possible insight into why the explicit schemes exhibit large amplitude instability.
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
179
The question remains as to whether the linear stability results can be interpreted in terms of the existence or otherwise of the energy conserving properties of scheme (A). A close investigation of the derivation of the stability results for scheme (A1) shows that the advantage arises from the term 1 + b2 τ 2 in the denominator of the coefficient of ξ in its characteristic equation. This term entered the equation because of the scheme’s implicit nature and is typical of the stability advantages conferred by implicit schemes. 4. A conservative explicit scheme and a non-conservative implicit scheme To gain evidence for the conjecture that it is the implicit nature of the schemes rather than an energy conservation property that defends against large amplitude instability we investigate a numerical scheme proposed by Fei and Vazquez [6] for the sine-Gordon equation. This scheme is easily generalized for the Klein–Gordon equation: (E) Generalization of the Fei and Vazquez sine-Gordon scheme [6]: n+1 F (un+2 1 1 2 n+2 j ) − F (uj ) n+3 n+2 n+1 n n+1 u − u + u + u − δ u + u + = 0. j j j j 2τ 2 j 2χ 2 x j un+2 − un+1 j j
The fact that this scheme is being used for an equation for which it was not strictly designed is not a drawback for our purposes. The scheme is second order accurate and, by Theorem 3, has a discrete conserved quantity associated with it. (E1) Application of the generalization of Fei and Vazquez scheme to the LKGE: 1 1 2 n+2 n+3 n+2 n+1 n n+1 2 n+2 n+1 u − u + u + u − δ u + u + b u + u =0 j j j j x j j j j 2τ 2 2χ 2 with von Neumann stability condition ρ 2 < 1 − 12 b2 τ 2 . This result shows us that the stability limit of the numerical scheme (E1) is dependent on the magnitude of b2 . In particular, for stability we require b2 < 2(1 − p 2 )/τ 2 . The scheme 1 2 n 1 2 n 1 n+1 n−1 (F) δ u − δ u + f u + f u =0 t j x j j j τ2 χ2 2 is implicit and has no known conserved energy. The von Neumann analysis of the associated linear problem gives results identical to those for scheme (A1). This suggests that this scheme will have a defense against large amplitude instability. The experiment carried out by Jiminez and Vazquez was repeated by applying the schemes (E) and (F) to the problem as described in their paper. (This experiment is described in Section 1.1 of this paper.) The nonlinear solve for scheme F was carried out using pointwise Newton–Raphson with iteration to convergence. The initial condition used in the experiments is repeated here for convenience: 2π U (0, x) = A 1 + cos x , 1.28 where A is a parameter. Scheme (F) showed no instability for this initial condition with the amplitude parameter A taking values up to 200, with runtime 200 and time step size 0.01. On the other hand, scheme (F) exhibits instability when the parameter A is as small as 24 (no observed instability when A = 23).
180
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
Table 1 A summary of the investigations carried out in this paper Scheme
LKGE stability limit
Exhibits large amplitude instability
Known conserved quantity
Implicit
(A)
ρ 2 < 1 + 14 b2 τ 2
No
Yes
Yes
(B)
ρ2
Yes
No
No
Yes
No
No
Yes
Yes
No
No
No
Yes
(C) and (D)
<1−
ρ2 < 1 − <1−
(E)
ρ2
(F)
ρ2 < 1 +
1 2 2 4b τ 1 2 2 4b τ 1 2 2 2b τ 1 2 2 4b τ
Both schemes were also tested with random initial conditions replacing the cosine term of the form = u1j (= u2j for scheme (E)) = A(2Uj ) where Uj is a random number in the range (0, 1). The results of this experiment confirmed the observations above. The solution to scheme (E) exhibited large amplitude instability for values of A as small as 20, whereas scheme (F) maintained its amplitude for values of A up to 300 (the limit of the test). Finally, the Jiminez and Vazquez experiment was repeated with A fixed at 60 and the A fixed at 120. In each case the step size was successively reduced by a half, starting at τ = 0.01, until no instability was observed during a run of 20,000 time steps. When A = 60 no instability was observed when τ = 0.0025 and when A = 120 no instability was observed when τ = 0.00125. This lends support to the view of the relationship between τ , A and instability suggested by the linear analysis carried out in Section 3, in the sense that large amplitude instability can defended against by a reduction in the time step size. u0j
5. Conclusion We summarize the results of our investigation in Table 1. The column titled ‘LKGE stability limit’ is the linear stability limit resulting from a von Neumann analysis on the scheme when it is applied to the linear Klein–Gordon equation (3.1)—see Section 3. The scheme (E) seems to break the link between large amplitude instability and the non-existence of a discrete conserved quantity for the numerical scheme. Scheme (F) gives further evidence to support the proposition that it is the implicit nature, rather than a conservation property, of the schemes that protects them from this type of instability. In this paper, the idea that the success of conservative finite difference schemes approximating the nonlinear Klein–Gordon equation can in part be measured by their defense against large amplitude instability has been questioned. The source of numerical instability for linear problems was used to gain insight into the behaviour of the nonlinear problems. From this analysis it has been argued that the defense against large amplitude instability displayed by nonlinear schemes is a consequence of their implicit nature.
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
181
Appendix For a list of definitions of the symbols used in this appendix see Section 2. Lemma 1. h1t uni , δt2 uni i = 1t /2h1t /2 uni , 1t /2 uni i. Proof.
1t unj , δt2 unj
n+1 = un+1 − un−1 − 2unj + un−1 j j , uj j
n+1 n n+1 n−1 n+1 n n−1 n−1 = un+1 − 2 un+1 − un−1 + 2 un−1 j , uj j , uj + uj , uj j , uj j , uj − uj , uj
(The middle two terms cancel. For ease of development a zero term is added to the end of the following equation.)
n+1/2
, uj
n+1 n n−1 n n−1 n−1 n n−1 n = un+1 − un+1 − un+1 j , uj j , uj + uj , uj − uj , uj j , uj + uj , uj
+ unj , unj − unj , unj = 1t /2 1t /2 unj , 1t /2unj .
2
Lemma 2. h1t unj , unj i = 1t /2 hunj , unj i. Proof.
X
n 1t unj , unj = un+1 − un−1 j j , uj = χ
n+1/2 n−1/2 uj
1t /2 uj
= 1t /2 uj
n−1/2
. 2
Lemma 3. hvjn , δx2 unj i = −h1Fx vjn , 1Fx unj i. Proof.
vjn , δx2 unj = vjn , unj+1 − 2unj + unj−1 = vjn , unj+1 − unj − vjn , unj − unj−1
= vjn , unj+1 − unj − vjn+1 , unj+1 − unj = − 1Fx vjn , 1Fx unj .
2
References [1] M.J. Ablowitz, M.D. Kruskal, J.F. Ladik, Solitary wave collisions, SIAM J. Appl. Math. 36 (1979) 428–437. [2] J. Argyris, M. Haase, J.C. Heinrich, Finite element approximations to the two-dimensional sine-Gordon solitons, Comput. Methods Appl. Mech. Engrg. 86 (1991) 1–26. [3] G. Ben-Yu, P.J. Pascual, M.J. Rodriguez, L. Vazquez, Numerical solution of the sine-Gordon equation, Appl. Math. Comput. 18 (1986) 1–14. [4] G. Ben-Yu, L. Xun, L. Vazquez, A Legendre spectral method for solving the nonlinear Klein–Gordon equation, Comput. Appl. Math. 15 (1) (1996) 19–36. [5] H. Elzoheiry, L. Iskandar, M.Sh. El-Deen Mohamedien, Iterative implicit schemes for the two and three dimensional sine-Gordon equation, J. Comput. Appl. Math. 34 (1991) 161–170. [6] Z. Fei, L. Vazquez, Two energy conserving numerical schemes for the sine-Gordon equation, Appl. Math. Comput. 35 (1990) 61–94. [7] S. Jiminez, Derivation of the discrete conservation laws of a family of finite difference schemes, Appl. Math. Comput. 64 (1994) 13–45.
182
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
[8] S. Jiminez, L. Vazquez, Analysis of four numerical schemes for a nonlinear Klein–Gordon equation, Appl. Math. Comput. 35 (1990) 61–94. [9] X. Li, B.Y. Guo, A Legendre spectral method for solving the nonlinear Klein–Gordon equation, J. Comput. Math. 15 (2) (1997) 105–126. [10] W.A. Strauss, L. Vazquez, Numerical solution of a nonlinear Klein–Gordon equation, J. Comput. Phys. 28 (1978) 271–278. [11] R. Vichnevetsky, J.B. Bowles, Fourier Analysis of Numerical Approximations to Hyperbolic Equations, SIAM Studies in Applied Mathematics (SIAM, Philadelphia, PA, 1982). [12] L. Vu-Quoc, S. Li, Invariant-conserving finite difference algorithms for the nonlinear Klein–Gordon equation, Comput. Methods Appl. Mech. Engrg. 107 (1993) 314–391. [13] Y.S. Wong, Q. Chang, L. Gong, An initial-boundary value problem of a nonlinear Klein–Gordon equation, Appl. Math. Comput. 84 (1997) 77–93.
Large amplitude instability in finite difference approximations to the Klein–Gordon equation Mark A.M. Lynch 1 Glasgow Caledonian University, Mathematics Department, City Campus, Cowcaddens Road, Glasgow G4 0BA, Scotland, UK
Abstract In a paper by Jiminez and Vazquez (1990), four finite difference schemes for approximating the nonlinear Klein– Gordon equation were discussed. They observed undesirable characteristics in some of the numerical schemes, in particular a loss of spatial symmetry and the onset of instability for large values of a parameter in the initial condition of the equation. The proposition, put forward by Jiminez and Vazquez, that the observed instability is protected against by using energy conserving schemes is questioned. An analysis of the schemes as applied to a linear problem is carried out and these indicate that the instability arises from the use of explicit finite difference schemes rather than any failure of energy conservation. This conjecture is further supported by an analysis of two further schemes. Some general results applicable to the nonlinear finite difference schemes are also presented. 1999 Elsevier Science B.V. and IMACS. All rights reserved. Keywords: Finite difference; Klein–Gordon; Stability
1. Introduction The numerical solutions to the NonLinear Klein–Gordon Equation (NLKGE) Ut t − Uxx + f (U ) = 0
(1.1)
have received considerable attention in the literature. A variety of second order finite difference schemes have been presented [8,12,13], and a strong emphasis has been given to developing schemes which conserve a discrete energy [7,8]. The latter approach has the particular appeal in that it reflects the energy conservation property of the differential equation. We note that alternative approaches using spectral and pseudospectral methods have recently been presented [4,9]. Investigations have been carried out on the particular case of the NLKGE, called the sine-Gordon equation [3,6], and, in particular, a finite element method is described in [2]. 1 E-mail: [email protected].
0168-9274/99/$20.00 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 8 ) 0 0 1 2 8 - 7
174
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
In a paper by Jiminez and Vazquez [8], four finite difference schemes for approximating the nonlinear Klein–Gordon equation were discussed. They observed undesirable characteristics in some of the numerical schemes, in particular a loss of spatial symmetry and the onset of instability for large values of a parameter in the equation’s initial condition. This instability will be referred to as large amplitude instability in this paper. The proposition, put forward by Jiminez and Vazquez [8], that the onset of this type of instability is protected against by using energy conserving schemes is questioned. In this paper the schemes are applied to the Linear Klein–Gordon Equation (LKGE) and a standard von Neumann analysis is performed to gain insight into the source of this large amplitude instability. The analysis suggests the counter proposition that the protection arises through the implicit nature of the schemes. To test this proposition two further schemes are considered: an explicit scheme with a known conserved energy and an implicit scheme with no known conserved energy. The explicit scheme used is an extension of the one described by Fei and Vazquez [6] for investigating the sine-Gordon equation. The von Neumann analysis suggests that the explicit scheme is not protected against the onset of the large amplitude instability, even though it has a conserved quantity, whereas the implicit scheme is protected from this type of instability. This conclusion is supported by numerical experiments with the two schemes. Before proceeding to an investigation of the particular schemes discussed by Jiminez and Vazquez, some general results are presented. These results cover the accuracy of some of the energy conserving schemes discussed in the literature and sufficiency conditions for the existence of an energy type conservation property. 1.1. Background We introduce some basic notation. Define τ as the time increment and χ is the spacing between points on the mesh (assumed uniform). Let U (t, x) be an exact solution to the NLKGE (1.1). Define a discrete function unj which is an approximation to the exact solution of the NLKGE, i.e., unj ≈ U (nτ, j χ). The schemes discussed in this paper will make use of the following central difference operators: δt2 unj = un+1 − 2unj + un−1 and δx2 unj = unj+1 − 2unj + unj−1 . Finally, we introduce the function F (U ), j j a primitive of the function f (U ). In [8] four second order finite difference (FD) schemes were subjected to a variety of numerical experiments. We concentrate on one of these experiments where the large amplitude of the initial condition resulted in three of the schemes failing to successfully propagate the solution whereas an energy conserving scheme succeeded. It is worth noting that none of the schemes exhibited any instability when the initial condition was of a ‘reasonable’ size. The schemes were set to solve the following initial value problem with periodic boundary condition: Ut t − Uxx + U + U 3 = 0 subject to the initial conditions
2π U (0, x) = A 1 + cos x 1.28
(1.2)
,
Ut (0, x) = 0,
defined on the space interval [0, 1.28]. The time step size was 0.01 and the mesh spacing 0.02. The parameter A was the subject of experimental variation. The four schemes used in the experimentation were:
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
175
(A) due to Strauss and Vazquez [10]: n−1 1 2 n 1 2 n F (un+1 j ) − F (uj ) δ u − δ u + = 0, τ2 t j χ2 x j un+1 − un−1 j j
(B) in common use: 1 2 n 1 δt uj − 2 δx2 unj + f unj = 0, 2 τ χ
(C) due to Ablowitz et al. [1]:
unj+1 + unj−1 1 2 n 1 2 n δ u − δ u +f τ2 t j χ2 x j 2
= 0,
(D) due to Jiminez [7]: 1 2 n 1 2 n F (unj+1 ) − F (unj−1 ) δ u − δ u + = 0. τ2 t j χ2 x j unj+1 − unj−1 The conclusions of the paper were that scheme (A) was superior to the other three in that it better preserved the spatial symmetry of the numerical solution and was not prone to the onset of instability for large values of the amplitude parameter A. Unfortunately, when looking for numerical schemes with higher orders of accuracy, it is not clear how to take advantage of this claim unless a conserved discrete energy can be found. In this paper we argue that this line of attack is not of itself a useful criterion for high order schemes. We use an analysis of the linear problem to propose that it is not the existence of the energy conservation property of the scheme that gives it protection from large amplitude instability, rather it is the implicit nature of the scheme.
2. Accuracy and a conserved quantity The energy conserving scheme discussed in [8] belongs to a family of schemes with the following structure: n+1/2
n−1/2
) − N(uj N(uj 1 2 n 1 δt uj − 2 δt2 unj + n+1 2 τ χ uj − un−1 j
)
= 0.
(2.1)
The symbol N represents some combination of the functions f and F with space and time differencing operators. We give two useful results relating to this generic scheme. Let u and v be discrete square summable functions defined for all integers (Z). Then hu, vi =
X
ui vi χ.
i∈Z
We define the following additional time differencing operators: 1t unj = un+1 − un−1 j j ,
n+1/2
Γt /2un = uj
and the following space differencing operators:
n−1/2
− uj
,
δt2 unj = un+1 − 2unj + un−1 j j ,
176
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
1Fx uj = unj+1 − unj , δx2 unj = unj+1 − 2unj + unj−1 , unj+1/2 + unj−1/2 µx uj = , 1x uj = unj+1 − unj−1 . 2 The application of space differencing operators to the functions u and v will be assumed to generate functions that are always square summable. Theorem 1. The numerical schemes (2.1) are second order accurate approximations to the nonlinear Klein–Gordon equation (1.1) if they satisfy n+1/2
n−1/2
) − N(uj ) ∂ (2.2) = 2 F unj + O τ 2 , τ ∂t where (∂/∂t)F (unj ) is interpreted as the time derivative of F evaluated at an exact solution of the partial differential equation at time j τ and at position j χ . N(uj
Proof. The linear part of scheme (2.1) is already a second order accurate approximation to the linear wave equation. Assume the scheme satisfies (2.2) this implies that n+1/2
n+1/2
) − N(uj ) ∂ n = 2f unj uj + O τ 2 . τ ∂t Using the second order central difference approximation to the time derivative of u and multiplying through by τ implies N(uj
n+1/2
N uj
n−1/2
− N uj
= un+1 − un−1 f unj + O τ 3 . j j
Dividing through by the order τ term un+1 − un−1 completes the proof. j j
2
Theorem 2. Solutions to the finite difference schemes (2.1) satisfy the following conservation condition:
n+1/2
1t /2 unj , 1t /2unj + ρ 2 1Fx uj
−1/2
, 1Fx uj
+ τ 2χ
X
N unj = K,
j
where K is a constant and the mesh ratio ρ is defined by ρ = τ/χ . Proof. Premultiply (2.1) by τ 2 χ1t unj and summing over all integer values j gives
n+1/2
1t unj , δt2 unj + 1t unj , ρ 2 δx2 unj + τ 2 1t unj ,
N(uj
n−1/2
) − N(uj
un+1 − un−1 j j
Application of the Lemmas 1–3 in Appendix reduces this to
1t /2
1t /2 unj , 1t /2 unj
+ρ
2
1Fx unj , 1Fx unj
+τ χ 2
X
N
unj
)
= 0.
= 0.
j
Hence the quantity in the square bracket is constant with respect to time. It is noted that for the particular NLKGE (1.2) the term N(u) will be positive definite and the conserved quantity is a true energy whenever ρ < 1 (see also the remarks by Jiminez [7]).
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182 n+1/2
177
n−1/2
Two particular choices of note are: N(unj ) = (F (uj ) + F (uj )) which returns scheme (A) n+1/2 n−1/2 n above, and N(uj ) = 2F ((uj + uj )/2) which returns a scheme described in [12] which is also a generalization of a scheme described in [6] for the sine-Gordon equation. 2 In Section 4 we will make use of the finite difference scheme described in Theorem 3. Theorem 3. Solutions to the finite difference schemes n+1 N(un+2 1 1 2 n+2 j ) − N(uj ) n+3 n+2 n+1 n n+1 u − uj + uj + uj − 2 δx uj + uj + =0 2τ 2 j 2χ un+2 − un+1 j j
satisfy the following conservation condition:
n+3/2
1t /2 uj
n+1/2
, 1t /2uj
n+3/2
+ ρ 2 1Fx uj
n+3/2
, 1Fx uj
+ 2τ 2 χ
X
n+3/2
N uj
=K
j
for some constant K. Proof. Multiply the difference equation through by 2τ 2 χ(un+2 − un+1 j j ), substitute ρ = τ/χ and sum over all integer values of j . Treating each of the three parts under the summation separately. The time differencing operator:
n+3 n+2 un+2 − un+1 − un+2 + un+1 + unj = un+2 − un+2 − un+1 j j , uj j j j j , 1t uj j
n+3/2
n+1/2
= 1t /2 1t /2 uj , 1t /2uj (applying Lemma 2 in Appendix). The space differencing operator:
2 n+2 un+2 − un+1 + δx2 un+1 j j , δx uj j
2 n+2 2 n+1 2 n+2 2 n+1 = un+2 − un+1 − un+1 + un+2 j , δx uj j , δx uj j , δx uj j , δx uj 2 n+2 2 n+1 = un+2 − un+1 j , δx uj j , δx uj
F n+2 F n+1 = − 1Fx un+2 + 1Fx un+1 j , 1x uj j , 1x uj
=
n+3/2 n+3/2 −1t /2 1Fx uj , 1Fx uj .
The nonlinear term is straightforward.
(applying Lemma 3 in Appendix)
2
It is an immediate consequence of Theorems 1 and 2 that the schemes (A)–(D) are second order and that scheme (A) has a conserved energy. The schemes (B)–(D) do not satisfy the conditions of Theorem 2 but this of itself does not mean that they fail to have a conserved quantity. As far as is known these schemes do not have a conserved energy, though scheme (D) (as does (A)) has a conserved momentum [8]. 3. Linear schemes The nonlinear schemes (A)–(D) are applied to the simpler linear Klein–Gordon equation Ut t − Uxx + b2 U = 0,
(3.1)
178
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
where b is a constant. The resulting schemes are then subjected to a von Neumann stability analysis. The resulting schemes are: (A1) (B1)
1 2 n 1 2 n b2 n+1 n−1 δ u − δ u + u + u = 0, t j x j j j τ2 χ2 2 1 2 n 1 δt uj − 2 δx2 unj + b2 unj = 0, 2 τ χ
(C1/D1) (note that (C) and (D) give rise to the same linear scheme): 1 2 n 1 2 n b2 n δ u − δx uj + uj +1 + unj−1 = 0. t j 2 2 τ χ 2
Applying a standard von Neumann stability analysis to the schemes (A1), (B1) and (C1/D1) results in characteristic equations with the following structure: ξ 2 − λξ + 1 = 0, where λ is determined by the scheme. The characteristic equations have two roots which we signify by ξ1 and ξ2 or simply ξ1,2 . Definition. The linear numerical schemes (A1), (B1), (C1/D1) are stable if and only if |ξ1,2 | 6 1 and ξ1 6= ξ2 whenever |ξ1,2 | = 1. The von Neumann analysis returns the following stability conditions on the mesh ratio ρ = τ/χ : • scheme (A1): ρ 2 < 1 + 14 b2 τ 2 , • scheme (B1) and (C1/D1): ρ 2 < 1 − 14 b2 τ 2 . The above results can be used to gain insight into the instability observed by Jiminez and Vazquez [8]. They observed in their numerical experiments that when the parameter A in the initial condition was allowed to become too large then the schemes (B)–(D) became unstable. Whereas scheme (A) remained stable for all observed values of the parameter A. An inspection of the above results shows that the stability of the schemes (B1) and (C1/D1) are adversely dependent on the parameter b2 , in the sense that increasing the value of b2 (with fixed τ ) eventually results in schemes that are unstable. On the other hand, the stability of the scheme (A1) is seen to be independent of the value of b2 . We note in passing that for both the schemes (B1) and (C1/D1) the onset of large amplitude instability can be compensated for by decreasing the time step size τ . The above arguments can be used to derive a possible explanation for the observations made by Jiminez and Vazquez in their numerical experiments. For large values of the parameter A in the initial condition the function f (U ) = U + U 3 will take values in the range 0 to order A3 . Any reasonable linear approximation of the form f (U ) ≈ b2 U will be a very steeply increasing function. The results of the preceding analysis indicate that when the schemes (B1) and (C1/D1) are applied to a linear function with sufficiently steep gradient numerical instability can be expected. On the other hand, no matter how large the value of b2 , the scheme (A1) will remain stable. This is not to say that other, nonlinear sources of instability, may be forced by the FD schemes but we have gained a possible insight into why the explicit schemes exhibit large amplitude instability.
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
179
The question remains as to whether the linear stability results can be interpreted in terms of the existence or otherwise of the energy conserving properties of scheme (A). A close investigation of the derivation of the stability results for scheme (A1) shows that the advantage arises from the term 1 + b2 τ 2 in the denominator of the coefficient of ξ in its characteristic equation. This term entered the equation because of the scheme’s implicit nature and is typical of the stability advantages conferred by implicit schemes. 4. A conservative explicit scheme and a non-conservative implicit scheme To gain evidence for the conjecture that it is the implicit nature of the schemes rather than an energy conservation property that defends against large amplitude instability we investigate a numerical scheme proposed by Fei and Vazquez [6] for the sine-Gordon equation. This scheme is easily generalized for the Klein–Gordon equation: (E) Generalization of the Fei and Vazquez sine-Gordon scheme [6]: n+1 F (un+2 1 1 2 n+2 j ) − F (uj ) n+3 n+2 n+1 n n+1 u − u + u + u − δ u + u + = 0. j j j j 2τ 2 j 2χ 2 x j un+2 − un+1 j j
The fact that this scheme is being used for an equation for which it was not strictly designed is not a drawback for our purposes. The scheme is second order accurate and, by Theorem 3, has a discrete conserved quantity associated with it. (E1) Application of the generalization of Fei and Vazquez scheme to the LKGE: 1 1 2 n+2 n+3 n+2 n+1 n n+1 2 n+2 n+1 u − u + u + u − δ u + u + b u + u =0 j j j j x j j j j 2τ 2 2χ 2 with von Neumann stability condition ρ 2 < 1 − 12 b2 τ 2 . This result shows us that the stability limit of the numerical scheme (E1) is dependent on the magnitude of b2 . In particular, for stability we require b2 < 2(1 − p 2 )/τ 2 . The scheme 1 2 n 1 2 n 1 n+1 n−1 (F) δ u − δ u + f u + f u =0 t j x j j j τ2 χ2 2 is implicit and has no known conserved energy. The von Neumann analysis of the associated linear problem gives results identical to those for scheme (A1). This suggests that this scheme will have a defense against large amplitude instability. The experiment carried out by Jiminez and Vazquez was repeated by applying the schemes (E) and (F) to the problem as described in their paper. (This experiment is described in Section 1.1 of this paper.) The nonlinear solve for scheme F was carried out using pointwise Newton–Raphson with iteration to convergence. The initial condition used in the experiments is repeated here for convenience: 2π U (0, x) = A 1 + cos x , 1.28 where A is a parameter. Scheme (F) showed no instability for this initial condition with the amplitude parameter A taking values up to 200, with runtime 200 and time step size 0.01. On the other hand, scheme (F) exhibits instability when the parameter A is as small as 24 (no observed instability when A = 23).
180
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
Table 1 A summary of the investigations carried out in this paper Scheme
LKGE stability limit
Exhibits large amplitude instability
Known conserved quantity
Implicit
(A)
ρ 2 < 1 + 14 b2 τ 2
No
Yes
Yes
(B)
ρ2
Yes
No
No
Yes
No
No
Yes
Yes
No
No
No
Yes
(C) and (D)
<1−
ρ2 < 1 − <1−
(E)
ρ2
(F)
ρ2 < 1 +
1 2 2 4b τ 1 2 2 4b τ 1 2 2 2b τ 1 2 2 4b τ
Both schemes were also tested with random initial conditions replacing the cosine term of the form = u1j (= u2j for scheme (E)) = A(2Uj ) where Uj is a random number in the range (0, 1). The results of this experiment confirmed the observations above. The solution to scheme (E) exhibited large amplitude instability for values of A as small as 20, whereas scheme (F) maintained its amplitude for values of A up to 300 (the limit of the test). Finally, the Jiminez and Vazquez experiment was repeated with A fixed at 60 and the A fixed at 120. In each case the step size was successively reduced by a half, starting at τ = 0.01, until no instability was observed during a run of 20,000 time steps. When A = 60 no instability was observed when τ = 0.0025 and when A = 120 no instability was observed when τ = 0.00125. This lends support to the view of the relationship between τ , A and instability suggested by the linear analysis carried out in Section 3, in the sense that large amplitude instability can defended against by a reduction in the time step size. u0j
5. Conclusion We summarize the results of our investigation in Table 1. The column titled ‘LKGE stability limit’ is the linear stability limit resulting from a von Neumann analysis on the scheme when it is applied to the linear Klein–Gordon equation (3.1)—see Section 3. The scheme (E) seems to break the link between large amplitude instability and the non-existence of a discrete conserved quantity for the numerical scheme. Scheme (F) gives further evidence to support the proposition that it is the implicit nature, rather than a conservation property, of the schemes that protects them from this type of instability. In this paper, the idea that the success of conservative finite difference schemes approximating the nonlinear Klein–Gordon equation can in part be measured by their defense against large amplitude instability has been questioned. The source of numerical instability for linear problems was used to gain insight into the behaviour of the nonlinear problems. From this analysis it has been argued that the defense against large amplitude instability displayed by nonlinear schemes is a consequence of their implicit nature.
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
181
Appendix For a list of definitions of the symbols used in this appendix see Section 2. Lemma 1. h1t uni , δt2 uni i = 1t /2h1t /2 uni , 1t /2 uni i. Proof.
1t unj , δt2 unj
n+1 = un+1 − un−1 − 2unj + un−1 j j , uj j
n+1 n n+1 n−1 n+1 n n−1 n−1 = un+1 − 2 un+1 − un−1 + 2 un−1 j , uj j , uj + uj , uj j , uj j , uj − uj , uj
(The middle two terms cancel. For ease of development a zero term is added to the end of the following equation.)
n+1/2
, uj
n+1 n n−1 n n−1 n−1 n n−1 n = un+1 − un+1 − un+1 j , uj j , uj + uj , uj − uj , uj j , uj + uj , uj
+ unj , unj − unj , unj = 1t /2 1t /2 unj , 1t /2unj .
2
Lemma 2. h1t unj , unj i = 1t /2 hunj , unj i. Proof.
X
n 1t unj , unj = un+1 − un−1 j j , uj = χ
n+1/2 n−1/2 uj
1t /2 uj
= 1t /2 uj
n−1/2
. 2
Lemma 3. hvjn , δx2 unj i = −h1Fx vjn , 1Fx unj i. Proof.
vjn , δx2 unj = vjn , unj+1 − 2unj + unj−1 = vjn , unj+1 − unj − vjn , unj − unj−1
= vjn , unj+1 − unj − vjn+1 , unj+1 − unj = − 1Fx vjn , 1Fx unj .
2
References [1] M.J. Ablowitz, M.D. Kruskal, J.F. Ladik, Solitary wave collisions, SIAM J. Appl. Math. 36 (1979) 428–437. [2] J. Argyris, M. Haase, J.C. Heinrich, Finite element approximations to the two-dimensional sine-Gordon solitons, Comput. Methods Appl. Mech. Engrg. 86 (1991) 1–26. [3] G. Ben-Yu, P.J. Pascual, M.J. Rodriguez, L. Vazquez, Numerical solution of the sine-Gordon equation, Appl. Math. Comput. 18 (1986) 1–14. [4] G. Ben-Yu, L. Xun, L. Vazquez, A Legendre spectral method for solving the nonlinear Klein–Gordon equation, Comput. Appl. Math. 15 (1) (1996) 19–36. [5] H. Elzoheiry, L. Iskandar, M.Sh. El-Deen Mohamedien, Iterative implicit schemes for the two and three dimensional sine-Gordon equation, J. Comput. Appl. Math. 34 (1991) 161–170. [6] Z. Fei, L. Vazquez, Two energy conserving numerical schemes for the sine-Gordon equation, Appl. Math. Comput. 35 (1990) 61–94. [7] S. Jiminez, Derivation of the discrete conservation laws of a family of finite difference schemes, Appl. Math. Comput. 64 (1994) 13–45.
182
M.A.M. Lynch / Applied Numerical Mathematics 31 (1999) 173–182
[8] S. Jiminez, L. Vazquez, Analysis of four numerical schemes for a nonlinear Klein–Gordon equation, Appl. Math. Comput. 35 (1990) 61–94. [9] X. Li, B.Y. Guo, A Legendre spectral method for solving the nonlinear Klein–Gordon equation, J. Comput. Math. 15 (2) (1997) 105–126. [10] W.A. Strauss, L. Vazquez, Numerical solution of a nonlinear Klein–Gordon equation, J. Comput. Phys. 28 (1978) 271–278. [11] R. Vichnevetsky, J.B. Bowles, Fourier Analysis of Numerical Approximations to Hyperbolic Equations, SIAM Studies in Applied Mathematics (SIAM, Philadelphia, PA, 1982). [12] L. Vu-Quoc, S. Li, Invariant-conserving finite difference algorithms for the nonlinear Klein–Gordon equation, Comput. Methods Appl. Mech. Engrg. 107 (1993) 314–391. [13] Y.S. Wong, Q. Chang, L. Gong, An initial-boundary value problem of a nonlinear Klein–Gordon equation, Appl. Math. Comput. 84 (1997) 77–93.