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

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.


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], selfinduced 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].
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 Ω = {(, ) ∈ [ 0  ,  1  ] × [ 0  ,  1  ]} for  > 0 and the parameter  being the so-called dissipative term, which is assumed to be a real number with  ≥ 0. When  = 0, (1) reduces to the undamped sine-Gordon equation in two space variables, whereas, when  > 0 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  2 / 2 , 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 (, )/ =  2 (, )/ 2 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 onedimensional 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.
Similarly, it is easily obtained that the truncation error of the CNFD scheme is O( 2 + ℎ 2  + ℎ 2  ).
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.

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 || ≤ 1 has to be satisfied.
The right-hand side of (20) gives if  = 0, then ( 23) is always satisfied, while when  > 0,  must satisfy the restriction condition Thus, the time step  needs to satisfy the stability conditions ( 22) and ( 24).

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: (34) The proof of Lemma 5 is complete.

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.

The Predictor of DFFD.
Using an analogous scheme as in [26,27], the predictor value û+1 , was evaluated from the following three-time level explicit scheme of DFFD ( 8): for ,  = 0, 1, . . ., .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.

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  = 0, (51) has the following form: (,  = 1, 2, . . ., ) ; (57) to approximate  −1 (, ) in the internal points the initial velocity is used, for this we discretize the initial velocity as Combining ( 57) with (58), we have , sin  0 , , ,  = 1, 2, . . .,  − 1.
In (64) and (65) the corrected values  +1 , instead of the predicted û+1 , 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.

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

Example 1.
To observe the behavior of the numerical method, let (, ) = 1 in (2); it is tested on the following problem: with initial conditions  (68) The theoretical solution of this problem, in which the parameter  = 0, is given by  (, , ) = 4 tan −1 (exp ( +  − )) . (69) As mentioned in Section 2, the truncation error of the DFFD scheme is ( 2 + ℎ 2  + ℎ 2  + (/ℎ  ) 2 + (/ℎ  ) 2 ).One can easily see that if we choose the time step  = ℎ 2  = ℎ 2  , for DFFD scheme and CNFD scheme, they are secondorder 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 ‖  −   ‖  2 , ‖  −   ‖  ∞ of the difference between the exact solution   given by (69) and the numerical solution   , the time is  = 7.0.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.
Take the same step length ℎ  = ℎ  = 0.25,  = 0.01 as in [21]; we apply the PC-DFFD scheme to computation of the solution to (66); the absolute errors are given at times  = 4, 7, 10, 20, 25. Figure 1 shows that when the space and time step are the same, compared with the absolute error at  = 7 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  = 25 or much longer time, where the absolute error is nearly the same as that at  = 4. 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 longtime numerical simulations.
It is known (see e.g., [7-9, 14-16, 18-21, 26, 27]) that, when  = 0, 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  = 0.001, ℎ  = ℎ  = 0.
When  = 0, 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  = 0 to  = 4, the solutions are almost not varying until while at  = 7 a deformation has appeared, while at  = 15 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  = 1.5, Figure 4 shows that at  = 9, even  = 15, the function graph virtually maintains the primary condition at  = 0, 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].

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 twodimensional sine-Gordon equation.Circular ring solitons are found for the case (, ) = 1 and initial conditions [12,14].Consider (73) In Figure 5, the numerical solutions of circular ring solitons for the  = 0 at  = 2.8, 5.6, 8.4, 11.2, 15, 18, 20 are shown in terms of sin(/2) for three-dimensional picture.The soliton from its initial position, where it appears as two homocentric ring solitons, is shrinking until  = 2.8 appears as a single-ring soliton.From  = 5.6, which could be considered as the beginning of the expansion phase, a radiation appears.This expansion continues until  = 11.2,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  = 0.05, 0.5, in Figure 6, with solution changes of a smaller  = 0.05 (Figure 6(a)) and a larger  = 0.5 (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.(

Example 4 (Collision of Two Circular Solitons
Numerical simulation presented in Figure 7 is for sin(/2) at levels  = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20 with  = 0, respectively.The solution which is shown includes the extension across  = −10 and  = −7 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

Figure 1 :
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  = 7 for test problems (66)-(68).

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

Table 3 :
The energy of the circular ring solitons.