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 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 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 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 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 . 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 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.
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.
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.
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].
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.
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).