Discrete Dynamics in Nature and Society
Volume 2017 (2017), Article ID 3634815, 8 pages
https://doi.org/10.1155/2017/3634815

Research Article

Compact Implicit Integration Factor Method for the Nonlinear Dirac Equation

Jing-Jing Zhang, Xiang-Gui Li, and Jing-Fang Shao

School of Applied Sciences, Beijing Information Science and Technology University, Beijing 100192, China

Correspondence should be addressed to Jing-Jing Zhang; jingjingzhang1112@163.com

Received 18 June 2017; Revised 7 September 2017; Accepted 19 September 2017; Published 18 October 2017

Academic Editor: Francisco R. Villatoro

Copyright © 2017 Jing-Jing Zhang 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 high-order accuracy numerical method is proposed to solve the -dimensional nonlinear Dirac equation in this work. We construct the compact finite difference scheme for the spatial discretization and obtain a nonlinear ordinary differential system. For the temporal discretization, the implicit integration factor method is applied to deal with the nonlinear system. We therefore develop two implicit integration factor numerical schemes with full discretization, one of which can achieve fourth-order accuracy in both space and time. Numerical results are given to validate the accuracy of these schemes and to study the interaction dynamics of the nonlinear Dirac solitary waves.

1. Introduction

It is well known that the Dirac equation plays a fundamental role in relativistic quantum physics [1]. As a relativistic wave equation, the Dirac equation can be used to describe the fine structure of the hydrogen atom and predict the existence of a positron verified by Anderson in 1933 [2]. The classical Dirac equation has been extensively known and studied since the graphene appeared in the laboratory in 2003 [3].

In the wake of developments in science and technology, the -dimensional nonlinear Dirac (NLD) equations have emerged as useful models in many physical problems, such as Bose-Einstein condensates in honeycomb optical lattices [4] and the gap solitons in nonlinear optics [5]. An important character of the NLD equation is that it has a solitary wave solution or particle-like solution. In order to construct the NLD model and keep its Lorentz invariance, the nonlinear self-interaction term can be built up in the Dirac Lagrangian function by using the bilinear covariants [6]. In addition, different self-interaction terms result in different NLD models, all of which have aroused extensive interest of many experts and scholars [7].

In the past few decades, many researchers have devoted themselves to studying the NLD equation analytically or numerically. Some analytical solitary wave solutions have been obtained in a series of research works [8–10]. With regard to more general initial conditions, however, it is difficult to derive the exact solutions of the NLD equation. Therefore, numerical techniques become important to solve the NLD equation and analyze its properties. From the numerical point of view, Alvarez et al. proposed a second-order Crank-Nicholson (CN) scheme [11] and the linearized CN scheme [12]. In 1989, split-step spectral schemes were presented by de Frutos and Sanz-Serna [13]. Modified Legendre rational spectral methods were developed for solving the -dimensional NLD equation by Wang and Guo in 2004 [14]. In 2006, Shao and Tang proposed a fourth-order accurate Runge-Kutta discontinuous Galerkin method [15]. In 2013, Xu et al. summarized several often-used finite deference schemes and applied the exponential operator splitting scheme to the NLD equation [7]. The CN finite difference method and the exponential wave integrator Fourier pseudospectral method were developed by Bao et al. in 2016 [16]. Besides, there are some other numerical techniques, such as multisymplectic Runge-Kutta methods [17], an efficient adaptive mesh redistribution method [18], an integrating factor method [19], and time-splitting methods with charge conservation [20]. Some experts also discussed the spin-orbit-coupled Bose-Einstein condensates [21] and the Maxwell-Dirac system [22–24] which are related to the NLD equation.

Recently, there has been a worldwide wave of interest in efficient and high-order numerical simulations for the NLD equation. Motivated by the previous work, we propose two compact implicit integration factor (IIF) schemes with high accuracy to solve the -dimensional NLD equation numerically. To this end, the fourth-order compact finite difference (CM) scheme is applied to discretize the NLD equation in spatial direction. Then, we adopt the classical IIF method [25] for the time evolution. A remarkable feature of the IIF method is that the nonlinear implicit terms are free of the exponential operator of the linear terms. Accordingly, when the IIF method is applied to the NLD equation, the size of the nonlinear system only relies on the number of the original equations. The method proposed can achieve fourth order in space and high order (arbitrary) in time [26, 27].

This paper is organized as follows. In Section 2, the -dimensional NLD equation and its two exact solutions are introduced in detail. In Section 3, we give the spatial semidiscrete scheme by using the fourth-order CM method. The IIF method is applied to the semidiscrete system in Section 4. In Section 5, numerical experiments are conducted to test the accuracy order of the schemes proposed. Moreover, the interaction dynamics of the solitary waves are discussed in this section. Conclusions are drawn in Section 6.

2. Nonlinear Dirac Equation and Its Two Exact Solutions

The NLD equation can be written in the covariant form as follows [7]:

where is the imaginary unit, is the rest mass, , is the transpose, and are complex-valued functions, represents the adjoint spinor, the superscript denotes the conjugate transpose, stands for the self-interaction Lagrangian, represents the covariant derivative , and denotes the Dirac matrices (); that is,

In this paper, we mainly concentrate our attention on the Soler model [28] with the following scalar self-interaction:

where is a real constant. For convenience, Pauli matrices are introduced in this paper. Let

and then substitute (3) into (1). We can derive the following -dimensional NLD equation:

The initial and Dirichlet boundary conditions of the NLD equation (5) are given as follows:

In addition, we can find that the -dimensional NLD equation (5) has two exact solutions [15, 29, 30], which are widely used in numerical experiments. The first is a standing wave solution which can be written as follows:

where

with .

The second is a single solitary wave solution which has the form

where denotes the wave velocity, , and is the sign function which is defined as

3. The Spatial Discretization

In this section, we start constructing the CM scheme which is an efficient and high-order accuracy method to discretize the -dimensional NLD equation (5) in the spatial direction.

In numerical computation, we cut the infinite interval into a finite interval which is proper and large enough. Thus, the original problem with the initial and Dirichlet boundary conditions can be approximated by

where is a positive constant.

Before giving the CM scheme, we first divide the domain uniformly with grid points and , where denotes the spatial step size, denotes the temporal step, and and are two positive integers.

Let and represent the numerical solution and the exact solution at , respectively. Next, we introduce some finite difference operators:

where can be obtained easily by Taylor series expansion [31]. Omitting the small terms , we can obtain the approximation of ; that is,

In the following, applying the above finite difference operator to problem (11), we can get the semidiscrete CM scheme which is a nonlinear ordinary differential system:

Hereafter, setting ,

we can rewrite the semidiscrete system (14) in the following matrix form:

4. The Temporal Discretization

To deal with the nonlinear system arising from the spatial discretization, we apply the IIF method which is an efficient time discretization method in this section. The method will be described in detail.

The operator on the right-hand side of the semidiscrete form (16) can be split into two parts, that is, and . To use the integration factor approach to (16), we multiply (16) by exponential matrix to derive

and then we integrate (17) over one time step from to to get

where .

Different approximation methods for the integral in (18) can produce different numerical methods. With the term unchanged, we can construct the second-order implicit exponential time differencing (ETD2) scheme by using the interpolation polynomial containing the points and to approximate . Denote as the numerical approximation of , and we can obtain the ETD2 scheme as follows:

Obviously, the system couples the linear term and the nonlinear term through . It is expensive to solve (11) by using the ETD2 scheme. In order to reduce the calculation cost greatly, we construct a class of IIF schemes which is decoupling between the implicit treatment of and the exact evaluation of . Setting

we will give a brief introduction to this IIF method.

A crucial step of constructing IIF schemes is to use Lagrange interpolation polynomial approximating the integral term in (18). To build an implicit scheme of the th order truncation error, we usually apply the th-order Lagrange interpolation polynomial including the point to approach the integral term . In this paper, we mainly restrict our attention to the construction of the second-order implicit scheme, CM-IIF2 scheme, and fourth-order implicit scheme, CM-IIF4 scheme. Hence, we set

Replacing the integral term in (18) with , we get the following CM-IIF2 scheme:

which is of fourth-order accuracy in space and of second-order accuracy in time. Similarly, to construct the CM-IIF4 scheme, we substitute for in (18) and obtain

The distinctive feature of the CM-IIF2 scheme and the CM-IIF4 scheme is that the exponential matrix is independent of and is only multiplied by known quantities. Therefore, the size of the nonlinear system is irrelevant to the number of spatial grid points, and it only hinges on the number of the original equations. In this paper, the term is firstly calculated with a given by using the function “expm” in Matlab and then stored for later use. However, for high-dimensional spatial discretization, the matrix becomes very large, and we can use the Krylov subspace approximations [26, 32] to the matrix exponential operator.

Because schemes (22)-(23) are implicit, an iterative algorithm is used to solve the problem in practical calculation. The iterative schemes can be rewritten as

In addition, the iteration termination condition for the above numerical techniques can be set as [33]

where is the number of iterations, represents a given error bound, and

5. Numerical Results

In this section, numerical examples are given to validate the accuracy of the numerical schemes proposed in Section 4, that is, CM-IIF2 scheme and CM-IIF4 scheme. And then the above schemes are applied to study the dynamic behaviors for the NLD solitary waves. All codes in this paper are run on an HP desktop computer with Intel Xeon CPU E5-1620 and 8 GB of RAM.

To quantify the numerical solution, the maximum error between the exact solution and the numerical solution at is defined as the discrete norm; that is,

5.1. Test of Accuracy

Example 1. We consider the following -dimensional NLD problem with zero Dirichlet boundary conditions:

where the initial data to this problem are derived by setting the rest mass in (9).

To test the accuracy of the CM-IIF2 scheme in space, we take the time size in numerical computation. Error analysis for the CM-IIF2 scheme at time is shown in Table 1. Similarly, we set the mesh size to validate the accuracy of the CM-IIF2 scheme in time. In Table 2, we give the maximum error between the exact solution and the numerical solution, error ratio, and accuracy order. As can be seen from Tables 1 and 2, the CM-IIF2 scheme can achieve fourth order in space and second order in time.

Table 1: The error, error ratio, and accuracy order computed by the CM-IIF2 scheme with at time (Example 1).

Error ratio
Order

0.4

—
—

0.2

16.5531
4.06

0.1

16.0127
4.00

0.05

16.3825
4.05

Table 2: The error, error ratio, and accuracy order computed by the CM-IIF2 scheme with at time (Example 1).

Error ratio
Order

0.1

—
—

0.05

4.0068
2.00

0.025

4.0096
2.00

0.0125

4.0324
2.01

In order to improve the computing efficiency, another method is used to test the accuracy of the CM-IIF4 scheme. Error analysis at time is shown in Table 3. We can draw a conclusion that the CM-IIF4 numerical scheme can achieve fourth order in both space and time.

Table 3: The error, error ratio, and accuracy order computed by the CT-IIF4 scheme at time (Example 1).

Mesh size

Error ratio
Order

—
—

18.0569
4.25

16.4403
4.05

15.7683
3.97

For the purpose of comparison, we give the following CN scheme which is of second-order accuracy in both space and time; that is,

where and the operators are defined by

The initial condition and boundary condition are defined as those of the above schemes. Using the CN scheme to calculate Example 1, we obtain some results which are listed in Table 4. From Table 4, we can see that the CN scheme achieves second order in both space and time. Therefore, the CM-IIF4 scheme proposed in this paper can not only improve the spatial accuracy but also improve the temporal accuracy.

Table 4: The error, error ratio, and accuracy order computed by the CN scheme at time (Example 1).

Mesh size

Error ratio
Order

—
—

3.8419
1.96

3.9990
2.00

3.9983
2.00

5.2. Interaction Dynamics

Example 2. The second example is to study the interaction dynamics of two moving solitary waves by solving (5). The initial data are given as follows:

where stands for the solitary wave placed at with parameters and represents the solitary wave placed at with .

In practical calculation, the problem is solved by using the CM-IIF2 scheme in a large domain . The time evolutions of the charge density from to are shown in Figure 1. The collision phenomenon of two solitary waves is depicted in Figure 1, from which we can see that the solitary waves keep moving with their original velocities and shapes after the collision [34–36].

Figure 1: The time evolution of the charge density at by using the CM-IIF2 scheme with and (Example 2).

Similarly, we further use the CM-IIF4 scheme to simulate the oscillating states. The result is shown in Figure 2. It can be seen that the elastic interaction happens for the two solitary waves.

Figure 2: The time evolution of the charge density at by using the CM-IIF4 scheme with and (Example 2).

6. Conclusion

In this paper, we propose a high-order accuracy numerical method for solving the -dimensional NLD equation. A compact spatial discretization method is proposed in this paper. Then, we apply the IIF method to treat the nonlinear system arising from the spatial discretization. In theory, the new compact IIF method can achieve fourth-order accuracy in space and arbitrary-order accuracy in time. In this paper, we mainly introduce the CM-IIF4 scheme which is of fourth order in both time and space. Numerical examples are devoted to confirming the view and investigating the dynamic behaviors of two NLD solitary waves.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was partially supported by grants from the National Natural Science Foundation of China (Project no. 11671044), the Science Challenge Project (Project no. TZ2016001), and the Beijing Municipal Education Commission (Project no. PXM2017_014224_000020).

References

P. A. M. Dirac, “The quantum theory of the electron,” Proceedings of the Royal Society of London A, vol. 117, no. 778, pp. 610–624, 1928.

C. D. Anderson, “The positive electron,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 43, no. 6, pp. 491–494, 1933.

K. S. Novoselov, A. K. Geim, S. V. Morozov et al., “Two-dimensional gas of massless Dirac fermions in graphene,” Nature, vol. 438, no. 7065, pp. 197–200, 2005.

L. H. Haddad and L. D. Carr, “The nonlinear Dirac equation in Bose-Einstein condensates: foundation and symmetries,” Physica D: Nonlinear Phenomena, vol. 238, no. 15, pp. 1413–1421, 2009.

I. V. Barashenkov, D. E. Pelinovsky, and E. V. Zemlyanaya, “Vibrations and oscillatory instabilities of gap solitons,” Physical Review Letters, vol. 80, no. 23, pp. 5117–5120, 1998.

J. Xu, S. Shao, H. Tang, and D. Wei, “Multi-hump solitary waves of a nonlinear Dirac equation,” Communications in Mathematical Sciences, vol. 13, no. 5, pp. 1219–1242, 2015.

J. Xu, S. Shao, and H. Tang, “Numerical methods for nonlinear Dirac equation,” Journal of Computational Physics, vol. 245, pp. 131–149, 2013.

S. Y. Lee, T. K. Kuo, and A. Gavrielides, “Exact localized solutions of two-dimensional field theories of massive fermions with Fermi interactions,” Physical Review D Particles & Fields, vol. 12, no. 8, pp. 2249–2253, 1975.

P. Mathieu, “Soliton solutions for Dirac equations with homogeneous non-linearity in (1+1) dimensions,” Journal of Physics A: Mathematical and General, vol. 18, no. 16, pp. L1061–L1066, 1985.

S. Shao, N. R. Quintero, F. G. Mertens, F. Cooper, A. Khare, and A. Saxena, “Stability of solitary waves in the nonlinear Dirac equation with arbitrary nonlinearity,” Physical Review E, vol. 90, no. 3, Article ID 032915, 2014.

A. Alvarez, P. Y. Kuo, and L. Vzquez, “The numerical study of a nonlinear one-dimensional Dirac equation,” Applied Mathematics and Computation, vol. 13, no. 1-2, pp. 1–15, 1983.

A. Alvarez, “Linearized CRAnk-NICholson scheme for nonlinear DIRac equations [CORrected title: Linearized CRAnk-NIColson scheme for nonlinear DIRac equations],” Journal of Computational Physics, vol. 99, no. 2, pp. 348–350, 1992.

J. de Frutos and J. M. Sanz-Serna, “Split-step spectral schemes for nonlinear Dirac systems,” Journal of Computational Physics, vol. 83, no. 2, pp. 407–423, 1989.

Z.-q. Wang and B.-y. Guo, “Modified Legendre rational spectral method for the whole line,” Journal of Computational Mathematics, vol. 22, no. 3, pp. 457–474, 2004.

S. Shao and H. Tang, “Higher-order accurate Runge-Kutta discontinuous Galerkin methods for a nonlinear Dirac model,” Discrete and Continuous Dynamical Systems - Series B, vol. 6, no. 3, pp. 623–640, 2006.

W. Bao, Y. Cai, X. Jia, and J. Yin, “Error estimates of numerical methods for the nonlinear Dirac equation in the nonrelativistic limit regime,” Science China Mathematics, vol. 59, no. 8, pp. 1461–1494, 2016.

J. Hong and C. Li, “Multi-symplectic Runge-Kutta methods for nonlinear Dirac equations,” Journal of Computational Physics, vol. 211, no. 2, pp. 448–472, 2006.

H. Wang and H. Tang, “An efficient adaptive mesh redistribution method for a non-linear Dirac equation,” Journal of Computational Physics, vol. 222, no. 1, pp. 176–193, 2007.

F. de la Hoz and F. Vadillo, “An integrating factor for nonlinear Dirac equations,” Computer Physics Communications, vol. 181, no. 7, pp. 1195–1203, 2010.

S. C. Li, X. G. Li, and F. Y. Shi, “Time-splitting methods with charge conservation for the nonlinear Dirac equation,” Numerical Methods for Partial Differential Equations, vol. 33, no. 5, pp. 1582–1602, 2017.

W. Bao and Y. Cai, “Ground states and dynamics of spin-orbit-coupled Bose-Einstein condensates,” SIAM Journal on Applied Mathematics, vol. 75, no. 2, pp. 492–517, 2015.

W. Bao and X.-G. Li, “An efficient and stable numerical method for the Maxwell-Dirac system,” Journal of Computational Physics, vol. 199, no. 2, pp. 663–687, 2004.

X.-G. Li, C. K. Chan, and Y. Hou, “A numerical method with particle conservation for the Maxwell-Dirac system,” Applied Mathematics and Computation, vol. 216, no. 4, pp. 1096–1108, 2010.

Z. Huang, S. Jin, P. A. Markowich, C. Sparber, and C. Zheng, “A time-splitting spectral scheme for the Maxwell-Dirac system,” Journal of Computational Physics, vol. 208, no. 2, pp. 761–789, 2005.

Q. Nie, Y.-T. Zhang, and R. Zhao, “Efficient semi-implicit schemes for stiff systems,” Journal of Computational Physics, vol. 214, no. 2, pp. 521–537, 2006.

S. Chen and Y.-T. Zhang, “Krylov implicit integration factor methods for spatial discretization on high dimensional unstructured meshes: application to discontinuous Galerkin methods,” Journal of Computational Physics, vol. 230, no. 11, pp. 4336–4352, 2011.

T. Jiang and Y.-T. Zhang, “Krylov implicit integration factor WENO methods for semilinear and fully nonlinear advection-diffusion-reaction equations,” Journal of Computational Physics, vol. 253, pp. 368–388, 2013.

M. Soler, “Classical, stable, Classical, stable, nonlinear spinor field with positive rest energy,” Physical Review D Particles & Fields, vol. 1, no. 10, pp. 2766–2769, 1970.

A. Alvarez and B. Carreras, “Interaction dynamics for the solitary waves of a nonlinear Dirac model,” Physics Letters A, vol. 86, no. 6-7, pp. 327–332, 1981.

S. Shao and H. Tang, “Interaction of solitary waves with a phase shift in a nonlinear Dirac model,” Communications in Computational Physics, vol. 3, no. 4, pp. 950–967, 2008.

T. Wang, B. Guo, and Q. Xu, “Fourth-order compact and energy conservative difference schemes for the nonlinear Schrödinger equation in two dimensions,” Journal of Computational Physics, vol. 243, pp. 382–399, 2013.

X.-G. Li, J. Zhu, R.-P. Zhang, and S. Cao, “A combined discontinuous Galerkin method for the dipolar Bose-Einstein condensation,” Journal of Computational Physics, vol. 275, pp. 363–376, 2014.

S. C. Li, X. G. Li, J. J. Cao, and W. B. Li, “High-order numerical method for the derivative nonlinear Schrödinger equation,” International Journal of Modeling, Simulation & Scientific Computing, vol. 8, Article ID 1750017, 2017.

G. Berkolaiko and A. Comech, “On spectral stability of solitary waves of nonlinear dirac equation in 1D,” Mathematical Modelling of Natural Phenomena, vol. 7, no. 2, pp. 13–31, 2012.

A. Contreras, D. E. Pelinovsky, and Y. Shimabukuro, “L2 orbital Stability of Dirac Solitons in the Massive Thirring Model,” Communications in Partial Differential Equations, vol. 41, no. 2, pp. 227–255, 2016.

J. Cuevas-Maraver, P. G. Kevrekidis, A. Saxena, A. Comech, and R. Lan, “Stability of solitary waves and vortices in a 2D nonlinear DIRac model,” Physical Review Letters, vol. 116, no. 21, Article ID 214101, 214101, 6 pages, 2016.