Discrete Dynamics in Nature and Society
Volume 2014 (2014), Article ID 784387, 22 pages
http://dx.doi.org/10.1155/2014/784387

Research Article

A Dufort-Frankel Difference Scheme for Two-Dimensional Sine-Gordon Equation

Zongqi Liang,¹ Yubin Yan,² and Guorong Cai¹

¹School of Sciences, Jimei University, Xiamen 361021, China
²Department of Mathematics, University of Chester, Chester CHI 4BJ, UK

Received 11 March 2014; Revised 22 July 2014; Accepted 24 July 2014; Published 29 October 2014

Academic Editor: Victor S. Kozyakin

Copyright © 2014 Zongqi Liang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A standard Crank-Nicolson finite-difference scheme and a Dufort-Frankel finite-difference scheme are introduced to solve two-dimensional damped and undamped sine-Gordon equations. The stability and convergence of the numerical methods are considered. To avoid solving the nonlinear system, the predictor-corrector techniques are applied in the numerical methods. Numerical examples are given to show that the numerical results are consistent with the theoretical results.

1. Introduction

The sine-Gordon equation arises in extended rectangular Josephson junctions, which consist of two layers of superconducting materials separated by an isolating barrier. A typical arrangement is a layer of lead and a layer of niobium separated by a layer of niobium oxide. A quantum particle has a nonzero significant probability of being able to penetrate in the other side of a potential barrier that would be impenetrable to the corresponding classical particle. Mathematical model of light bullets in Maxwell-Bloch system can also be described by the two-dimensional undamped sine-Gordon equation [1]. This equation is also applied in a large number of areas of physics, for example, crystal dislocation theory [2], self-induced transparency [2], laser physics [2], and particle physics [3, 4]. It is also a special case of the Bady Skyrme model which describes baryons in a nonlinear manner [5].

Like some well-known partial differential equations such as KdV, mKdV, and nonlinear Schrödinger equations [6], the two-dimensional sine-Gordon equation possesses various types of soliton solutions (such as line solitons, elliptic ring soliton, and ring solitons). For the undamped () sine-Gordon equation in higher dimensions exact soliton solutions have been obtained in Hirota [7], in Zagrodzinsky [8] using Lamb’s method, in Leibbrandt [9] by Bläcklund transformation, and in Kaliappan and Lakshmanan [10] by Painlevé transcendents.

Numerical solutions for two-dimensional undamped sine-Gordon equation have been given among others by Guo et al. [11] using two finite difference schemes, Xin [1] studying sine-Gordon equation as an asymptotic reduction of the two level dissipationless Maxwell-Bloch system, Christiansen and Lomdahl [12] using a generalized leapfrog method, Argyris et al. [13] using the finite element method. Sheng et al. [14] introduced a split cosine scheme, Bratsos [15] used a three-time level fourth-order explicit finite-difference scheme, Mirzaei and Dehghan [16] applied the continuous linear boundary elements method, Chen et al. [17] applied the multilevel augmentation method for solving the sine-Gordon equations, and so forth. Numerical approaches to the damped sine-Gordon equation can also be found in Nakajima et al. [18] who considered dimensionless loss factors and unitless normalized bias and Gorria et al. [19] who studied the nonlinear wave propagation in a planar wave guide consisting of two rectangular regions joined by a bent of constant curvature. Bratsos [15] and Djidjeli et al. [20] used a two-step one-parameter leap-frog scheme, which is a generalization to that used by Christiansen and Lomdahl [12]. Dehghan and Shokri [21] used the radial basis functions as a truly meshfree method, to solve the two-dimensional damped/undamped sine-Gordon equation.

Consider the two-dimensional sine-Gordon equation as follows:

with in the region for and the parameter being the so-called dissipative term, which is assumed to be a real number with . When , (1) reduces to the undamped sine-Gordon equation in two space variables, whereas, when to the damped one. The function can be explained as a Josephson current density.

Initial conditions associated with (1) will be assumed to be of the form

with initial velocity

In (2) and (3), the functions and represent wave modes or kinks and velocity, respectively. Boundary conditions will be assumed to be of the form

where and are normal gradients along the boundary of the region .

If partial differential equation contains the second order term , the equation is discredited by using the standard Crank-Nicolson scheme. However, it is verified that the scheme is unstable unconditionally for the heat equation when it is discretized by the leap-frog scheme in time direction and the Crank-Nicolson scheme in space direction, see [22]. The Dufort-Frankel method is very similar to the leap-frog scheme but it has better numerical stability, and sometimes it will bring unpredicted effect (if Dufort-Frankel approximate scheme is applied in spatial direction, the heat equation is unconditionally stable [22]). Wu [23] applied the Dufort-Frankel scheme for solving the linear and nonlinear one-dimensional Schrödinger equations and obtained an unconditionally stable scheme. Markowich et al. [24] applied the Wigner-measure analysis for investigating the convergence of the Dufort-Frankel scheme for the Schrödinger equation in semiclassical regime. Lai et al. [25] used a simple Dufort-Frankel type scheme for solving the time-dependent Gross-Pitaevskii equation. In this paper, we will use the Dufort-Frankel scheme to solve two-dimensional sine-Gordon equation.

The organization of the paper is as follows. In Section 2, a Crank-Nicolson finite-difference (CNFD) scheme and a Dufort-Frankel finite-difference (DFFD) scheme for solving two-dimensional sine-Gordon equation are introduced. In Section 3, the stabilities of the two schemes are discussed. In Section 4, the error estimates of CNFD and DFFD schemes are proved. To avoid solving nonlinear equations, the predictor-corrector methods of the two schemes are proposed in Section 5. In Section 6, numerical results are investigated.

2. Numerical Methods

2.1. The CNFD Scheme and the DFFD Scheme

For the numerical solution the region with its boundary consisting of the lines and being covered with a rectangular mesh, the points with coordinates , with , and . The and represent the space steps along direction and direction, while represents the time step. The solution of an approximating difference scheme at the same point will be denoted by () for the purpose of analyzing stability, the numerical value actually obtained will be denoted by .

The CNFD scheme of (1) is

where , .

The DFFD scheme for (5) is obtained by replacing the term by in the CNFD scheme in the discretizations of and , that is,

Note that

where . The DFFD scheme can be expressed as follows, with :

2.2. Local Truncation Error

By Tayor expansion, it can be easily obtained that for the smooth function , the principal part of the local truncation error of the DFFD (8) is

To make (8) be consistent, it is known from (9) that the step sizes have to be limited such that when , , and , we have , . Letting , the above requirement holds automatically.

Similarly, it is easily obtained that the truncation error of the CNFD scheme is .

Remark 1. In terms of the compatibility conditions, the constraint of time and space steps in the DFFD scheme is more demanding than that of the CNFD scheme, but it is easy for the time step to satisfy it, which is a usual method. In the following discussion of convergence, it is found that the time and space steps in the CNFD scheme have to satisfy the constraint; however, the DFFD scheme is unconditional.

Remark 2. The DFFD is a scheme in which a term is added to the leading coefficient of DNFD. This extra term enhances the stability of it and greatly reduces the restriction on time step.

3. Stability Analysis of CNFD and DFFD

In this section, the stabilities of the two schemes are discussed. Following the Fourier method of analyzing stability a small error of the following form is considered:

with

where is a complex number and are real. Due to von Neumann criterion for stability, the condition has to be satisfied.

3.1. Stability Analysis for the DFFD

Using (10)-(11) and Maclaurin’s expansion

(9) can be written as

where

Here we linearize the term in square brackets and . Using (14) and , we get the following stability equation:

where

Assume that , if , , , then holds for any . Hence holds. Through von Neumann criterion the DFFD scheme is unconditionally stable.

Assume that , , and , and . Note that always holds. Hence the solutions,

of (15) satisfy , which implies that if , , the scheme is stable.

3.2. Stability Analysis for the CNFD

Similarly, the stability equation of CNFD (5) is

where

To guarantee , we need or

The left-hand side of (20) has the form

which is true if

Let be the roots of (22) and . Note that ; the roots are real and distinct. Let the positive root ; we then have .

The right-hand side of (20) gives

if , then (23) is always satisfied, while when , must satisfy the restriction condition

Thus, the time step needs to satisfy the stability conditions (22) and (24).

4. Convergence and Error Estimates for DFFD and CNFD

We define a discrete inner product and its associated norm by

Conveniently, we note that

The difference scheme DFFD (8) can be written in the following form:

We have the following Lemmas. Consider

Lemma 3. Consider

Lemma 4. Consider

Lemma 5. Letting , for arbitrary , we have, for ,

Proof. By inequality,

the definition of the difference quotient implies

By using the definition of norm, we get

Formulas (31)-(32) imply that Lemma 5 holds. In particular,

The proof of Lemma 5 is complete.

Next we consider the convergence of the DFFD (8). Suppose , we define the following weak form, for :

Subtracting (8) from (27), we have

where

In (36), letting ,

Applying Lemmas 3 and 4, we have

Summing up for from 1 to and multiplying , we get

where

Note that

which follows from

Similarly

From (42) and (44), we obtain

By Lemma 5, we have

in the last formula, we need to choose suitable to satisfy . If the time step satisfies , then we will make all coefficients bigger than 0 and obtain

By Gronwall’s inequality, we have, for all ,

Remark 6. The convergence analyses for the CNFD scheme are the same as for the DFFD scheme. The difference is that there is no in the first term of CNFD scheme (5); so in (46), the convergence of the CNFD scheme needs is satisfying both and for . We choose , the convergence of CDFD needs satisfying the constraint condition

5. The Predictor-Corrector Scheme

To avoid solving the nonlinear system arising from systems (5) and (8), the following predictor-corrector (P-C) scheme is used.

5.1. The Predictor of DFFD

Using an analogous scheme as in [26, 27], the predictor value was evaluated from the following three-time level explicit scheme of DFFD (8):

for .

Following a similar approach for the stability analysis of the nonlinear scheme as in Section 3.1, it can be proved that the characteristic equation of the predictor is given by

And the scheme (51) is therefore unconditionally stable.

5.2. The Predictor of CNFD

The predictor of CNFD satisfies

Following the discussion in Section 3.2, the characteristic equation of (53) is

Denote , the roots of (54) are

Let with . To ensure , we require , that is,

5.3. The Predictor Algorithm Implementation for the DFFD Scheme

The predictor of the DFFD scheme is a three-time level explicit scheme; in order to obtain the same second order accuracy, we deal with the initial value of the DFFD scheme, that is, at , (51) has the following form:

to approximate in the internal points the initial velocity is used, for this we discretize the initial velocity as

Combining (57) with (58), we have

Thus, (51) holds; applying the boundary conditions for and ,

The second order discretization form of the boundary condition is

Substituting (61) into (60), we have

Similarly, for and ,

5.4. The Correctors of DFFD and CNFD

The corrector of DFFD can be proposed as

for .

The corrector of CNFD can be proposed as

for .

In (64) and (65) the corrected values instead of the predicted were also used, and the stability analysis of the corrector is analogous to that developed in Sections 3.1 and 3.2. The initial value problem and boundary problem are solved in the same way as in Section 5.3.

6. Numerical Results

In this section we present some numerical results of two schemes for the two-dimensional sine-Gordon equation.

6.1. Example 1

To observe the behavior of the numerical method, let in (2); it is tested on the following problem:

with initial conditions

and boundary conditions

The theoretical solution of this problem, in which the parameter , is given by

As mentioned in Section 2, the truncation error of the DFFD scheme is . One can easily see that if we choose the time step , for DFFD scheme and CNFD scheme, they are second-order convergence. Respectively by the PC-DFFD scheme and the PC-CNFD scheme we take the same space step and work out equation (66), where the initial conditions (67) and the boundary conditions (68) are employed. The error is measured by the , of the difference between the exact solution given by (69) and the numerical solution , the time is . Table 1 shows that it needs a sufficiently small space and time step to keep the stability of PC-CNFD scheme while the PC-DFFD scheme is unconditionally stable, which is in consistent with the theoretical results.

Table 1: Two different errors for and at .

Take the same step length , as in [21]; we apply the PC-DFFD scheme to computation of the solution to (66); the absolute errors are given at times . Figure 1 shows that when the space and time step are the same, compared with the absolute error at in [21], the accuracy of the PC-DFFD scheme is much better than the method in [21] and the absolute error becomes about 0.1 times smaller than the RBF method in [21]. Compared with [13, 15, 18–21, 26, 27], the PC-DFFD scheme is much better than the numerical algorithms presented in these articles. Figure 2 shows that the PC-DFFD scheme has better stability at or much longer time, where the absolute error is nearly the same as that at . Compared with the scheme proposed in [15, 18–21, 26, 27], the PC-DFFD scheme maintains its simplicity and better stabilities. The accumulation of absolute errors can not lead to infinite increases of them; hence, the scheme can be applied in long-time numerical simulations.

(a)

(b)

Figure 1: (a) Absolute error of the solution of RBF method [21] and (b) absolute error of the solution of the PC-DFFD scheme in for test problems (66)–(68).

(a)

(b)

(c)

(d)

Figure 2: Absolute error of the solution of the PC-DFFD scheme in for test problems (66)–(68).

It is known (see e.g., [7–9, 14–16, 18–21, 26, 27]) that, when , for the sine-Gordon equation the energy given by the following expression,

is conserved. We also investigate this property for sine-Gordon equation; the evaluation of is performed using the composite trapezoidal rule for integration.

In the numerical calculations that follow, various cases involving line and ring solitons for the solution of (5) are reported. In all the following experiments, by the PC-DFFD scheme,we choose , , and the boundary conditions are taken to be

6.2. Example 2 (Superposition of Two Line Solitons)

It has been known that superpositions of two orthogonal line solitons can be acquired for , with initial conditions [14]

When , the results in Figure 3 show the break up of two orthogonal line solitons which are parallel to the diagonal and are moving away from each other in the direction of , undisturbed. From to , the solutions are almost not varying until while at a deformation has appeared, while at its deformation and original morphology will have a great change, but its energy remains basically the same. Table 2 shows that the scheme keeps better energy conservation and reliability, compared with others presented by Dehghan and Shokri [21] and Mirzaei and Dehghan [16], Bratsos [15, 26, 27]. For a small value of , the dissipative term is found to have little effect on the superposition of two line solitons. For a large value of , Figure 4 shows that at , even , the function graph virtually maintains the primary condition at , while its energy decreases as increases and damages its attributes of energy conservation illustrated in Table 2. However, the dissipative term is found to slow down the separation and break-up of two orthogonal line solitons as time increases. The present solutions are in good agreement with the corresponding results of [15, 16, 18–21, 26, 27].

Table 2: The energy of the superposition of two orthogonal line solitions.

Figure 3: The numerical solutions at times for superposition of two line solitons.

(a)

(b)

Figure 4: The numerical solutions at times , with for superposition of two line solitons.

6.3. Example 3 (Circular ring Solitons)

Bogolyubskii and Makhankov [28], Bogolyubskii [29], and Christiansen and Lomdahl [12] have investigated numerically the behavior of a circular ring quasisoliton or pulsion arising from the two-dimensional sine-Gordon equation. Circular ring solitons are found for the case and initial conditions [12, 14]. Consider

In Figure 5, the numerical solutions of circular ring solitons for the at are shown in terms of for three-dimensional picture. The soliton from its initial position, where it appears as two homocentric ring solitons, is shrinking until appears as a single-ring soliton. From , which could be considered as the beginning of the expansion phase, a radiation appears. This expansion continues until , where the soliton is almost reformed. These results are in agreement with the published ones in [15, 16, 18–21, 26, 27]. Table 3 presents the values of at some selected times , that shows the energy remains constant for a longer time as time increases. For the different , in Figure 6, with solution changes of a smaller (Figure 6(a)) and a larger (Figure 6(b)), the results are the same as [21]. As increases, the initial shrunk ring soliton was found to be changing more slowly from its initial position as time increases; the dissipative term is slowing down the evolution of the line soliton as time increases.

Table 3: The energy of the circular ring solitons.

Time

Figure 5: Circular ring solitons solution at with .

Figure 6: Circular ring solitons solution at with (a) and (b) .

6.4. Example 4 (Collision of Two Circular Solitons)

The collision of two expanding circular ring solitons is considered with and initial conditions

Numerical simulation presented in Figure 7 is for at levels = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20 with , respectively. The solution which is shown includes the extension across and by symmetry properties of the problem [14, 20, 21]. Figure 7 demonstrates the collision between two expanding circular ring solitons in which two smaller ring solitons bounding an annular region emerge into a large ring soliton. The simulated solution is again precisely consistent to existing results; contour maps are given to show more clearly the movement of solitons. Though minor disturbances can be observed in middle of the numerical solution, probably due to the transactions following the symmetry features in computations, the overall simulation results match well those described in [12, 14, 16, 20, 21] with satisfaction. Simultaneously, it is found that after , the originally formed large ring soliton is to split into two line solitons to the boundary extension.

Figure 7: Collision of two circular solitons, at with .

6.5. Example 5 (Collision of Four Circular ring Solitons)

A collision of four expanding circular ring solitons is investigated, for , and initial conditions [14]

The simulation is based on an extension across and due to the symmetry. The results are depicted in Figure 8 for at = 0, 2.5, 4, 5, 6, 7.5, 8, 10, 12, and 15 in terms of ; Figure 8 demonstrates precisely the collision between four expanding circular ring solitons in which the smaller ring solitons bounding an annular region emerge into a large ring soliton. The results are in good agreement with corresponding surfaces given in [13, 14, 16, 20, 21, 26]. Similarly, as time increases, at the large soliton is to split and form 5 or even more solutions of varying sizes.

Figure 8: Collision of four circular ring solitons, at with .

7. Conclusion

In this paper, two numerical methods, CNFD and DFFD, are constructed for solving the two-dimensional sin-Gordon equation. The stability, convergence, and error estimates are discussed. The DFFD scheme is unconditionally stable and convergent, while the CNFD scheme requires more critical space and time step constraints, in order to guarantee its stability and convergence. We establish two explicit PC-CNFD schemes and PC-DFFD scheme, whose stabilities are discussed. The PC-DFFD scheme is unconditionally stable. The numerical experiments indicate that the PC-DFFD scheme has better stability in comparison with the methods in [13, 15, 18–21, 26, 27] and is proper to the long-time computation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The research is supported by the Natural Science Foundation of Fujian Province (Grant no. 2012J01013, 2013J01245), the Higher College Special Foundation for Education Department of Fujian Province (Grant no. JK2012025), the Natural Science Foundation of Xiamen City (Grant no. 3502Z20110010), and the National Natural Science of China (Grant no. 11201178). The authors would like to thank Professor Chuanju Xu for insightful discussions regarding the two-dimensional sine-Gordon equation.

References

J. X. Xin, “Modeling light bullets with the two-dimensional sine-Gordon equation,” Physica D, vol. 135, no. 3-4, pp. 345–368, 2000.

L. Bergé, “Wave collapse in physics: principles and applications to light and plasma waves,” Physics Reports, vol. 303, no. 5-6, pp. 259–370, 1998.

L. Bergé and T. Colin, “A singular perturbation problem for an envelope equation in plasma physics,” Physica D, vol. 84, pp. 437–459, 1995.

N. Bloembergen, M. Mlejnek, J. Moloney, and E. Wright, “Femtosecond pulses at the boundary of a nonlinear dielectric,” in Advances in Laser Physics, V. S. Letokhov and P. Meystre, Eds., Gordon & Breach, New York, NY, USA, 1999.

T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Physical Review Letters, vol. 78, no. 17, pp. 3282–3286, 1997.

W. E. Fitzgibbon, “Strongly damped quasilinear evolution equations,” Journal of Mathematical Analysis and Applications, vol. 79, no. 2, pp. 536–550, 1981.

R. Hirota, “Exact three-soliton solution of the two-dimensional sine-Gordon equation,” Journal of the Physical Society of Japan, vol. 35, no. 5, article 1566, 1973.

J. Zagrodzinsky, “Particular solutions of the sine-Gordon equation in 2 + 1 dimensions,” Physics Letters A, vol. 72, no. 4-5, pp. 284–286, 1979.

G. Leibbrandt, “New exact solutions of the classical sine-Gordon equation in 2 + 1 and 3 + 1 dimensions,” Physical Review Letters, vol. 41, no. 7, pp. 435–438, 1978.

P. Kaliappan and M. Lakshmanan, “Kadomtsev-Petviashvili and two-dimensional sine-Gordon equations: reduction to Painlevé transcendents,” Journal of Physics A: Mathematical and General, vol. 12, pp. 233–249, 1979.

B. Y. Guo, P. J. Pascual, M. J. Rodriguez, and L. Vázquez, “Numerical solution of the sine-Gordon equation,” Applied Mathematics and Computation, vol. 18, no. 1, pp. 1–14, 1986.

P. L. Christiansen and P. S. Lomdahl, “Numerical study of 2 + 1 dimensional sine-Gordon solitons,” Physica D, vol. 2, no. 3, pp. 482–494, 1981.

J. Argyris, M. Haase, and J. C. Heinrich, “Finite element approximation to two-dimensional sine-Gordon solitons,” Computer Methods in Applied Mechanics and Engineering, vol. 86, no. 1, pp. 1–26, 1991.

Q. Sheng, A. Q. M. Khaliq, and D. A. Voss, “Numerical simulation of two-dimensional sine-Gordon solitons via a split cosine scheme,” Mathematics and Computers in Simulation, vol. 68, no. 4, pp. 355–373, 2005.

A. G. Bratsos, “A modified predictor-corrector scheme for the two-dimensional sine-Gordon equation,” Numerical Algorithms, vol. 43, no. 4, pp. 295–308, 2006.

D. Mirzaei and M. Dehghan, “Boundary element solution of the two-dimensional sine-Gordon equation using continuous linear elements,” Engineering Analysis with Boundary Elements, vol. 33, no. 1, pp. 12–24, 2009.

J. Chen, Z. Chen, and S. Cheng, “Multilevel augmentation methods for solving the sine-Gordon equation,” Journal of Mathematical Analysis and Applications, vol. 375, no. 2, pp. 706–724, 2011.

K. Nakajima, Y. Onodera, T. Nakamura, and R. Sato, “Numerical analysis of vortex motion on Josephson structures,” Journal of Applied Physics, vol. 45, no. 9, pp. 4095–4099, 1974.

C. Gorria, Y. B. Gaididei, M. P. Soerensen, P. L. Christiansen, and J. G. Caputo, “Kink propagation and trapping in a two-dimensional curved Josephson junction,” Physical Review B, vol. 69, Article ID 134506, 2004.

K. Djidjeli, W. G. Price, and E. H. Twizell, “Numerical solutions of a damped sine-Gordon equation in two space variables,” Journal of Engineering Mathematics, vol. 29, no. 4, pp. 347–369, 1995.

M. Dehghan and A. Shokri, “A numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 700–715, 2008.

K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equation, Cambridge University Press, Cambridge, UK, 2005.

L. Wu, “Dufort-Frankel-type methods for linear and nonlinear Schrödinger equations,” SIAM Journal on Numerical Analysis, vol. 33, no. 4, pp. 1526–1533, 1996.

P. A. Markowich, P. Pietra, and C. Pohl, “A Wigner-measure analysis of the Dufort-Frankel scheme for the Schrödinger equation,” SIAM Journal on Numerical Analysis, vol. 40, no. 4, pp. 1281–1310, 2002.

M.-C. Lai, C.-Y. Huang, and T.-S. Lin, “A simple Dufort-Frankel-type scheme for the Gross-Pitaevskii equation of Bose-Einstein condensates on different geometries,” Numerical Methods for Partial Differential Equations, vol. 20, no. 4, pp. 624–638, 2004.

A. G. Bratsos, “The solution of the two-dimensional sine-Gordon equation using the method of lines,” Journal of Computational and Applied Mathematics, vol. 206, no. 1, pp. 251–277, 2007.

A. G. Bratsos, “A third order numerical scheme for the two-dimensional sine-Gordon equation,” Mathematics and Computers in Simulation, vol. 76, no. 4, pp. 271–282, 2007.

I. L. Bogolyubskii and V. G. Makhankov, “Life time of pulsating solitons in certain classical models,” JETP Letters, vol. 24, no. 1, pp. 12–14, 1976.

I. L. Bogolyubskii, “Oscillating particle-like solutions of the nonlinear Klein-Gordon equation,” JETP Letters, vol. 24, no. 10, pp. 535–538, 1976.