Research Article  Open Access
JingJing Zhang, XiangGui Li, JingFang Shao, "Compact Implicit Integration Factor Method for the Nonlinear Dirac Equation", Discrete Dynamics in Nature and Society, vol. 2017, Article ID 3634815, 8 pages, 2017. https://doi.org/10.1155/2017/3634815
Compact Implicit Integration Factor Method for the Nonlinear Dirac Equation
Abstract
A highorder 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 fourthorder 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 BoseEinstein 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 particlelike solution. In order to construct the NLD model and keep its Lorentz invariance, the nonlinear selfinteraction term can be built up in the Dirac Lagrangian function by using the bilinear covariants [6]. In addition, different selfinteraction 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 secondorder CrankNicholson (CN) scheme [11] and the linearized CN scheme [12]. In 1989, splitstep spectral schemes were presented by de Frutos and SanzSerna [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 fourthorder accurate RungeKutta discontinuous Galerkin method [15]. In 2013, Xu et al. summarized several oftenused 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 RungeKutta methods [17], an efficient adaptive mesh redistribution method [18], an integrating factor method [19], and timesplitting methods with charge conservation [20]. Some experts also discussed the spinorbitcoupled BoseEinstein condensates [21] and the MaxwellDirac system [22–24] which are related to the NLD equation.
Recently, there has been a worldwide wave of interest in efficient and highorder 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 fourthorder 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 fourthorder 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 complexvalued functions, represents the adjoint spinor, the superscript denotes the conjugate transpose, stands for the selfinteraction 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 selfinteraction: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 highorder 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 bywhere 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 righthand 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 deriveand then we integrate (17) over one time step from to to getwhere .
Different approximation methods for the integral in (18) can produce different numerical methods. With the term unchanged, we can construct the secondorder 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 . Settingwe 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 thorder Lagrange interpolation polynomial including the point to approach the integral term . In this paper, we mainly restrict our attention to the construction of the secondorder implicit scheme, CMIIF2 scheme, and fourthorder implicit scheme, CMIIF4 scheme. Hence, we set
Replacing the integral term in (18) with , we get the following CMIIF2 scheme:which is of fourthorder accuracy in space and of secondorder accuracy in time. Similarly, to construct the CMIIF4 scheme, we substitute for in (18) and obtain
The distinctive feature of the CMIIF2 scheme and the CMIIF4 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 highdimensional 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, CMIIF2 scheme and CMIIF4 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 E51620 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 CMIIF2 scheme in space, we take the time size in numerical computation. Error analysis for the CMIIF2 scheme at time is shown in Table 1. Similarly, we set the mesh size to validate the accuracy of the CMIIF2 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 CMIIF2 scheme can achieve fourth order in space and second order in time.


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

For the purpose of comparison, we give the following CN scheme which is of secondorder 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 CMIIF4 scheme proposed in this paper can not only improve the spatial accuracy but also improve the temporal accuracy.

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 CMIIF2 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].
Similarly, we further use the CMIIF4 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.
6. Conclusion
In this paper, we propose a highorder 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 fourthorder accuracy in space and arbitraryorder accuracy in time. In this paper, we mainly introduce the CMIIF4 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. View at: Google Scholar
 C. D. Anderson, “The positive electron,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 43, no. 6, pp. 491–494, 1933. View at: Publisher Site  Google Scholar
 K. S. Novoselov, A. K. Geim, S. V. Morozov et al., “Twodimensional gas of massless Dirac fermions in graphene,” Nature, vol. 438, no. 7065, pp. 197–200, 2005. View at: Publisher Site  Google Scholar
 L. H. Haddad and L. D. Carr, “The nonlinear Dirac equation in BoseEinstein condensates: foundation and symmetries,” Physica D: Nonlinear Phenomena, vol. 238, no. 15, pp. 1413–1421, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar
 J. Xu, S. Shao, H. Tang, and D. Wei, “Multihump solitary waves of a nonlinear Dirac equation,” Communications in Mathematical Sciences, vol. 13, no. 5, pp. 1219–1242, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 J. Xu, S. Shao, and H. Tang, “Numerical methods for nonlinear Dirac equation,” Journal of Computational Physics, vol. 245, pp. 131–149, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 S. Y. Lee, T. K. Kuo, and A. Gavrielides, “Exact localized solutions of twodimensional field theories of massive fermions with Fermi interactions,” Physical Review D Particles & Fields, vol. 12, no. 8, pp. 2249–2253, 1975. View at: Google Scholar
 P. Mathieu, “Soliton solutions for Dirac equations with homogeneous nonlinearity in (1+1) dimensions,” Journal of Physics A: Mathematical and General, vol. 18, no. 16, pp. L1061–L1066, 1985. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Google Scholar
 A. Alvarez, P. Y. Kuo, and L. Vzquez, “The numerical study of a nonlinear onedimensional Dirac equation,” Applied Mathematics and Computation, vol. 13, no. 12, pp. 1–15, 1983. View at: Publisher Site  Google Scholar  MathSciNet
 A. Alvarez, “Linearized CRAnkNICholson scheme for nonlinear DIRac equations [CORrected title: Linearized CRAnkNIColson scheme for nonlinear DIRac equations],” Journal of Computational Physics, vol. 99, no. 2, pp. 348–350, 1992. View at: Publisher Site  Google Scholar  MathSciNet
 J. de Frutos and J. M. SanzSerna, “Splitstep spectral schemes for nonlinear Dirac systems,” Journal of Computational Physics, vol. 83, no. 2, pp. 407–423, 1989. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Google Scholar  MathSciNet
 S. Shao and H. Tang, “Higherorder accurate RungeKutta discontinuous Galerkin methods for a nonlinear Dirac model,” Discrete and Continuous Dynamical Systems  Series B, vol. 6, no. 3, pp. 623–640, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar  MathSciNet
 J. Hong and C. Li, “Multisymplectic RungeKutta methods for nonlinear Dirac equations,” Journal of Computational Physics, vol. 211, no. 2, pp. 448–472, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 H. Wang and H. Tang, “An efficient adaptive mesh redistribution method for a nonlinear Dirac equation,” Journal of Computational Physics, vol. 222, no. 1, pp. 176–193, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar  MathSciNet
 S. C. Li, X. G. Li, and F. Y. Shi, “Timesplitting methods with charge conservation for the nonlinear Dirac equation,” Numerical Methods for Partial Differential Equations, vol. 33, no. 5, pp. 1582–1602, 2017. View at: Google Scholar
 W. Bao and Y. Cai, “Ground states and dynamics of spinorbitcoupled BoseEinstein condensates,” SIAM Journal on Applied Mathematics, vol. 75, no. 2, pp. 492–517, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 W. Bao and X.G. Li, “An efficient and stable numerical method for the MaxwellDirac system,” Journal of Computational Physics, vol. 199, no. 2, pp. 663–687, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 X.G. Li, C. K. Chan, and Y. Hou, “A numerical method with particle conservation for the MaxwellDirac system,” Applied Mathematics and Computation, vol. 216, no. 4, pp. 1096–1108, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 Z. Huang, S. Jin, P. A. Markowich, C. Sparber, and C. Zheng, “A timesplitting spectral scheme for the MaxwellDirac system,” Journal of Computational Physics, vol. 208, no. 2, pp. 761–789, 2005. View at: Publisher Site  Google Scholar  MathSciNet
 Q. Nie, Y.T. Zhang, and R. Zhao, “Efficient semiimplicit schemes for stiff systems,” Journal of Computational Physics, vol. 214, no. 2, pp. 521–537, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar  MathSciNet
 T. Jiang and Y.T. Zhang, “Krylov implicit integration factor WENO methods for semilinear and fully nonlinear advectiondiffusionreaction equations,” Journal of Computational Physics, vol. 253, pp. 368–388, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Google Scholar
 A. Alvarez and B. Carreras, “Interaction dynamics for the solitary waves of a nonlinear Dirac model,” Physics Letters A, vol. 86, no. 67, pp. 327–332, 1981. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Google Scholar  MathSciNet
 T. Wang, B. Guo, and Q. Xu, “Fourthorder 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. View at: Publisher Site  Google Scholar  MathSciNet
 X.G. Li, J. Zhu, R.P. Zhang, and S. Cao, “A combined discontinuous Galerkin method for the dipolar BoseEinstein condensation,” Journal of Computational Physics, vol. 275, pp. 363–376, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 S. C. Li, X. G. Li, J. J. Cao, and W. B. Li, “Highorder numerical method for the derivative nonlinear Schrödinger equation,” International Journal of Modeling, Simulation & Scientific Computing, vol. 8, Article ID 1750017, 2017. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar  MathSciNet
 J. CuevasMaraver, 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. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2017 JingJing 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.