# Compact Implicit Integration Factor Method for the Nonlinear Dirac Equation

Authors: Jing-Jing Zhang; Xiang-Gui Li; Jing-Fang Shao
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A high-order accuracy numerical method is proposed to solve the(1+1)-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.


---

## Body

## 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(1+1)-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 (1+1)-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(1+1)-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 Section2, the (1+1)-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]:(1)iγμ∂μ-mΦ+∂L1Φ∂Φ¯=0,where i=-1 is the imaginary unit, m≥0 is the rest mass, Φ=Φ(x,t)=(ϕ1(x,t),ϕ2(x,t))†, † is the transpose, ϕ1 and ϕ2 are complex-valued functions, Φ¯:=Φ∗γ0 represents the adjoint spinor, the superscript ∗ denotes the conjugate transpose, L1[Φ] stands for the self-interaction Lagrangian, ∂μ represents the covariant derivative ∂0=∂t,∂1=∂x, and γμ denotes the Dirac matrices (μ=0,1); that is, (2)γ0=100-1,γ1=01-10.In this paper, we mainly concentrate our attention on the Soler model [28] with the following scalar self-interaction:(3)L1Φ=λΦ¯Φ2,where λ is a real constant. For convenience, Pauli matrices are introduced in this paper. Let (4)σ1=γ0γ1=0110,σ2=-iγ1,σ3=γ0, and then substitute (3) into (1). We can derive the following (1+1)-dimensional NLD equation:(5)iσ3∂tΦ-σ2∂xΦ-mΦ+2λΦ∗σ3ΦΦ=0,x∈R,t>0.The initial and Dirichlet boundary conditions of the NLD equation (5) are given as follows:(6)Φx,0=Φ0x=ϕ10x,ϕ20x†,x∈R,limΦx,tx→∞=0,x∈R,t>0.In addition, we can find that the(1+1)-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:(7)Φx,tsw=ϕ1swx,tϕ2swx,t=AxiBxe-iwt,where (8)Ax=1/λm2-w2m+wcoshm2-w2xm+wcosh2m2-w2x,Bx=1/λm2-w2m-wsinhm2-w2xm+wcosh2m2-w2x, with 0<w≤m.The second is a single solitary wave solution which has the form(9)ϕ1ssx,t=γ+12ϕ1swx~,t~+signvγ-12ϕ2swx~,t~,ϕ2ssx,t=signvγ-12ϕ1swx~,t~+γ+12ϕ2swx~,t~, where v denotes the wave velocity, γ=1/1-v2,x~=γ(x-vt),t~=γ(t-vx), and sign(x) is the sign function which is defined as (10)signx=-1,x<0,0,x=0,1,x>0.


## 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(1+1)-dimensional NLD equation (5) in the spatial direction.In numerical computation, we cut the infinite interval into a finite intervalΩ=(L1,L2)∈R which is proper and large enough. Thus, the original problem with the initial and Dirichlet boundary conditions can be approximated by(11)∂tΦ+σ1∂xΦ+imσ3Φ-2iλΦ∗σ3Φσ3Φ=0,x∈Ω,t∈0,T,Φx,0=Φ0x,x∈Ω¯,ΦL1,t=ΦL2,t=0,t∈0,T,where T is a positive constant.Before giving the CM scheme, we first divide the domain[L1,L2]×[0,T] uniformly with grid points xj=L1+jh(j=0,1,2,…,J+1) and tn=nτ(n=0,1,2,…,N), where h=L2-L1/J+1 denotes the spatial step size, τ=T/N denotes the temporal step, and J and N are two positive integers.LetΦj and Φ(xj) represent the numerical solution and the exact solution at xj, respectively. Next, we introduce some finite difference operators: (12)δxΦj=Φj+1-Φj-12h,ThΦj=16Φj+1+4Φj+Φj-1, where Th(∂Φ/∂x|(xj))=Φ(xj+1)-Φ(xj-1)/2h+O(h4) can be obtained easily by Taylor series expansion [31]. Omitting the small terms O(h4), we can obtain the approximation of ∂/∂xΦ(xj); that is, (13)∂∂xΦxj≈Th-1Φxj+1-Φxj-12h=Th-1δxΦj.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:(14)dΦjdt+Th-1σ1δxΦj+imσ3Φj-2iλΦj∗σ3Φjσ3Φj=0,j=1,2,…,J,Φ0=ΦJ+1=0.Hereafter, settingY=Φ1,Φ2,…,ΦJ†=ϕ11,ϕ21,ϕ12,ϕ22,…,ϕ1J,ϕ2J†,L=A-1Γ, (15)O2=0000,I2=1001,A=164I2I2I24I2I2⋱⋱⋱⋱⋱I2I24I2J×J,Γ=12hO2-σ1σ1O2-σ1⋱⋱⋱⋱⋱-σ1σ1O2J×J,NY=-mσ3+2λΦ1∗σ3Φ1σ3-mσ3+2λΦ2∗σ3Φ2σ3⋱-mσ3+2λΦJ∗σ3ΦJσ3J×J,we can rewrite the semidiscrete system (14) in the following matrix form:(16)dYdt=LY+iNYY.


## 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, LY and iN(Y)Y. To use the integration factor approach to (16), we multiply (16) by exponential matrix e-Lt to derive(17)e-LtdYdt=e-LtLY+ie-LtNYY,and then we integrate (17) over one time step from tn to tn+1=tn+τ to get(18)Ytn+1=eLτYtn+ieLτ∫0τe-LsNYtn+sYtn+sds,where 0≤s≤τ.Different approximation methods for the integral in (18) can produce different numerical methods. With the term e-Ls unchanged, we can construct the second-order implicit exponential time differencing (ETD2) scheme by using the interpolation polynomial containing the points tn+1 and tn to approximate NYtn+sYtn+s. Denote Yn as the numerical approximation of Ytn, and we can obtain the ETD2 scheme as follows: (19)Yn+1=eLτYn+ieLτ-I2LNYn+1Yn+1+NYnYn.Obviously, the system couples the linear term and the nonlinear term througheLτ. 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 N(Y)Y and the exact evaluation of LY. Setting(20)ps=e-LsNYtn+sYtn+s,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 termp(s) in (18). To build an implicit scheme of the rth order truncation error, we usually apply the r-1th-order Lagrange interpolation polynomial including the point tn+1 to approach the integral term p(s). 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(21)q1s=1τNYnYnτ-s+e-LτNYn+1Yn+1s,q2s=16τ3-e2LτNYn-2Yn-2s-τss+τ+3eLτNYn-1Yn-1s-τss+2τ-3NYnYns-τs+τs+2τ+e-LτNYn+1Yn+1ss+τs+2τ.Replacing the integral termps in (18) with q1s, we get the following CM-IIF2 scheme:(22)Yn+1=eLτYn+iτ2NYnYn+iτ2NYn+1Yn+1,which is of fourth-order accuracy in space and of second-order accuracy in time. Similarly, to construct the CM-IIF4 scheme, we substitute q2s for ps in (18) and obtain(23)Yn+1=eLτYn+1924iτNYnYn-524ieLττNYn-1Yn-1+124ie2LττNYn-2Yn-2+924iτNYn+1Yn+1.The distinctive feature of the CM-IIF2 scheme and the CM-IIF4 scheme is that the exponential matrixeLτ is independent of N(Yn+1) 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 eLτ 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 L 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 (24)Yn+1k+1=eLτYn+iτ2NYnYn+iτ2NYn+1kYn+1k,Yn+1k+1=eLτYn+1924iτNYnYn-524ieLττNYn-1Yn-1+124ie2LττNYn-2Yn-2+924iτNYn+1kYn+1k.In addition, the iteration termination condition for the above numerical techniques can be set as [33] (25)Yn+1k+1-Yn+1k<ϵ, where k is the number of iterations, ϵ=10-13 represents a given error bound, and (26)Yn+1k=0=Yn,limk→+∞Yn+1k+1=Yn+1,0≤n≤N.


## 5. Numerical Results

In this section, numerical examples are given to validate the accuracy of the numerical schemes proposed in Section4, 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Φ(xj,tn) and the numerical solution Φjn at t=tn is defined as the discrete L∞ norm; that is, (27)e∞tn≔max0≤j≤JΦxj,tn-Φjn.


### 5.1. Test of Accuracy

Example 1.
We consider the following(1+1)-dimensional NLD problem with zero Dirichlet boundary conditions:(28)∂tΦ+σ1∂xΦ+imσ3Φ-2iλΦ∗σ3Φσ3Φ=0,x∈-50,50,t>0,Φ-50,t=Φ50,t=0,t≥0,where the initial data to this problem are derived by setting the rest mass m=1,λ=1,w=4/5,v=1/10 in (9).To test the accuracy of the CM-IIF2 scheme in space, we take the time sizeτ=10-5 in numerical computation. Error analysis for the CM-IIF2 scheme at time t=1 is shown in Table 1. Similarly, we set the mesh size h=0.1 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τ=10-5 at time t=1 (Example 1).
h
e
∞
Error ratio
Order
0.4
4.9206
e
-
05
—
—
0.2
2.9726
e
-
06
16.5531
4.06
0.1
1.8564
e
-
07
16.0127
4.00
0.05
1.1332
e
-
08
16.3825
4.05Table 2
The error, error ratio, and accuracy order computed by the CM-IIF2 scheme withh=0.1 at time t=1 (Example 1).
τ
e
∞
Error ratio
Order
0.1
1.8281
e
-
04
—
—
0.05
4.5625
e
-
05
4.0068
2.00
0.025
1.1379
e
-
05
4.0096
2.00
0.0125
2.8219
e
-
06
4.0324
2.01In order to improve the computing efficiency, another method is used to test the accuracy of the CM-IIF4 scheme. Error analysis at timet=1 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 timet=1 (Example 1).
Mesh size
e
∞
Error ratio
Order
h
=
0.8
,
τ
=
0.04
8.9423
e
-
04
—
—
h
/
2
,
τ
/
2
4.9523
e
-
05
18.0569
4.25
h
/
4
,
τ
/
4
3.0123
e
-
06
16.4403
4.05
h
/
8
,
τ
/
8
1.9103
e
-
07
15.7683
3.97For the purpose of comparison, we give the following CN scheme which is of second-order accuracy in both space and time; that is,(29)δt+Φjn+σ1δxΦjn+1/2+imσ3Φjn+1/2-iλΦjn+1∗σ3Φjn+1+Φjn∗σ3Φjnσ3Φjn+1/2=0, where xj,tn∈[L1,L2]×[0,T] and the operators are defined by (30)Φjn+1/2=Φjn+1+Φjn2,δt+Φjn=Φjn+1-Φjnτ.The initial condition and boundary condition are defined as those of the above schemes. Using the CN scheme to calculate Example1, 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 timet=1 (Example 1).
Mesh size
e
∞
Error ratio
Order
h
=
0.2
,
τ
=
0.1
0.0011
—
—
h
/
2
,
τ
/
2
2.8632
e
-
04
3.8419
1.96
h
/
4
,
τ
/
4
7.1597
e
-
05
3.9990
2.00
h
/
8
,
τ
/
8
1.7907
e
-
05
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:(31)Φx,0=Φleftssx+6,0+Φrightssx-6,0,where Φleftssx+6,0 stands for the solitary wave placed at x=-6 with parameters v=4/7,w=4/5,m=λ=1 and Φrightssx-6,0 represents the solitary wave placed at x=6 with v=-9/10,w=4/5,m=λ=1.In practical calculation, the problem is solved by using the CM-IIF2 scheme in a large domain-50,50. The time evolutions of the charge density from t=0 to t=30 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 atx=0 by using the CM-IIF2 scheme with h=0.05 and τ=0.01 (Example 2).Similarly, we further use the CM-IIF4 scheme to simulate the oscillating states. The result is shown in Figure2. It can be seen that the elastic interaction happens for the two solitary waves.Figure 2
The time evolution of the charge density atx=0 by using the CM-IIF4 scheme with h=0.05 and τ=0.01 (Example 2).


## 5.1. Test of Accuracy

Example 1.
We consider the following(1+1)-dimensional NLD problem with zero Dirichlet boundary conditions:(28)∂tΦ+σ1∂xΦ+imσ3Φ-2iλΦ∗σ3Φσ3Φ=0,x∈-50,50,t>0,Φ-50,t=Φ50,t=0,t≥0,where the initial data to this problem are derived by setting the rest mass m=1,λ=1,w=4/5,v=1/10 in (9).To test the accuracy of the CM-IIF2 scheme in space, we take the time sizeτ=10-5 in numerical computation. Error analysis for the CM-IIF2 scheme at time t=1 is shown in Table 1. Similarly, we set the mesh size h=0.1 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τ=10-5 at time t=1 (Example 1).
h
e
∞
Error ratio
Order
0.4
4.9206
e
-
05
—
—
0.2
2.9726
e
-
06
16.5531
4.06
0.1
1.8564
e
-
07
16.0127
4.00
0.05
1.1332
e
-
08
16.3825
4.05Table 2
The error, error ratio, and accuracy order computed by the CM-IIF2 scheme withh=0.1 at time t=1 (Example 1).
τ
e
∞
Error ratio
Order
0.1
1.8281
e
-
04
—
—
0.05
4.5625
e
-
05
4.0068
2.00
0.025
1.1379
e
-
05
4.0096
2.00
0.0125
2.8219
e
-
06
4.0324
2.01In order to improve the computing efficiency, another method is used to test the accuracy of the CM-IIF4 scheme. Error analysis at timet=1 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 timet=1 (Example 1).
Mesh size
e
∞
Error ratio
Order
h
=
0.8
,
τ
=
0.04
8.9423
e
-
04
—
—
h
/
2
,
τ
/
2
4.9523
e
-
05
18.0569
4.25
h
/
4
,
τ
/
4
3.0123
e
-
06
16.4403
4.05
h
/
8
,
τ
/
8
1.9103
e
-
07
15.7683
3.97For the purpose of comparison, we give the following CN scheme which is of second-order accuracy in both space and time; that is,(29)δt+Φjn+σ1δxΦjn+1/2+imσ3Φjn+1/2-iλΦjn+1∗σ3Φjn+1+Φjn∗σ3Φjnσ3Φjn+1/2=0, where xj,tn∈[L1,L2]×[0,T] and the operators are defined by (30)Φjn+1/2=Φjn+1+Φjn2,δt+Φjn=Φjn+1-Φjnτ.The initial condition and boundary condition are defined as those of the above schemes. Using the CN scheme to calculate Example1, 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 timet=1 (Example 1).
Mesh size
e
∞
Error ratio
Order
h
=
0.2
,
τ
=
0.1
0.0011
—
—
h
/
2
,
τ
/
2
2.8632
e
-
04
3.8419
1.96
h
/
4
,
τ
/
4
7.1597
e
-
05
3.9990
2.00
h
/
8
,
τ
/
8
1.7907
e
-
05
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:(31)Φx,0=Φleftssx+6,0+Φrightssx-6,0,where Φleftssx+6,0 stands for the solitary wave placed at x=-6 with parameters v=4/7,w=4/5,m=λ=1 and Φrightssx-6,0 represents the solitary wave placed at x=6 with v=-9/10,w=4/5,m=λ=1.In practical calculation, the problem is solved by using the CM-IIF2 scheme in a large domain-50,50. The time evolutions of the charge density from t=0 to t=30 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 atx=0 by using the CM-IIF2 scheme with h=0.05 and τ=0.01 (Example 2).Similarly, we further use the CM-IIF4 scheme to simulate the oscillating states. The result is shown in Figure2. It can be seen that the elastic interaction happens for the two solitary waves.Figure 2
The time evolution of the charge density atx=0 by using the CM-IIF4 scheme with h=0.05 and τ=0.01 (Example 2).


## 6. Conclusion

In this paper, we propose a high-order accuracy numerical method for solving the(1+1)-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.



---

*Source: 3634815-2017-10-18.xml*