Mesoscopic Simulation of the Two-Component System of Coupled Sine-Gordon Equations with Lattice Boltzmann Method

In this paper, a new lattice Boltzmann model for the two-component system of coupled sine-Gordon equations is presented by using the coupled mesoscopic Boltzmann equations. Via the Chapman-Enskog multiscale expansion, the macroscopical governing evolution system can be recovered correctly by selecting suitable discrete equilibrium distribution functions and the amending functions. The mesoscopic model has been validated by several related issues where analytic solutions are available. The experimental results show that the numerical results are consistent with the analytic solutions. From the mesoscopic point of view, the present approach provides a new way for studying the complex nonlinear partial differential equations arising in natural nonlinear phenomena of engineering and science.


Introduction
It is well known that most of the nonlinear phenomena that arise in engineering fields and mathematical physics, including plasma physics, fluid dynamics and nonlinear fiber optics, can be described by nonlinear partial differential equations (NPDEs). NPDEs have become an available tool for describing these natural nonlinear phenomena of engineering and science models. Hence, it becomes more and more important to be acquainted with all traditional and recently developed methods for NPDEs, and implementation of these methods [1,2]. Some of the most interesting features or physical rules are concealed in their nonlinear characteristics and can only be researched with an approximate method that is designed for inherent nonlinearity issues. As a result of the complexity and nonlinearity of the wave evolution equations, there is no uniform approach to obtain all solutions of the nonlinear wave evolution system. Hence, to find more precise and more effective methods for acquiring the nonlinear wave evolution equations has been an attractive research business. In the last few decades, quite a number of research work has been designed to research various types of nonlinear wave evolution equations. They include effective and broadly applicable techniques such as the finite difference method, variational iteration method, finite element method, finite volume method, boundary elements method, etc.
In this work, we consider the two-component system of coupled sine-Gordon equations, which was introduced recently by Khusnutdinova and Pelinovsky [48]. The basic one-dimensional form is shown as follows: where α remarks the ratio of the acoustic velocities between the components u and w, the dimensionless parameter δ 2 is the same with the ratio of masses of particles in the "lower" and the "upper" parts of the crystal. This system produces the Frenkel-Kontorova dislocation model [49], and this system has also turned out to be highly suitable to describe fluxon phenomena of stacked intrinsic Josephson junctions in high temperature superconductors. Moreover, this system has been studied extensively for two-junction stacks, for stacks consisting of more junctions only some special cases have been analyzed [50]. In addition, this system (1) with α = 1 was proposed to describe the open states in DNA [51]. We consider the above system (1) with the initial conditions as follows: and the boundary conditions: where ϕ 1 (x), ψ 1 (x), ϕ 2 (x), ψ 2 (x) and φ i (t)(i = 1, 2, 3, 4) are known functions.
There are many analytical methods solving the two-component system of coupled sine-Gordon equations, such as the modified decomposition method [52], the homotopy analysis method [53], the hyperbolic auxiliary function method [54], the homotopy perturbation method [55], the rational exponential ansatz method [56], the variational iteration method [57] and the modified Kudryashov method [58]. However, to our best knowledge, there are few numerical method to solve this coupled system. In recent years, the studies in Refs. [33,59,60] show that the LB method may be an valid numerical solver for real and complex nonlinear coupled systems. Therefore, it is worthy to more study LB method and enlarge its applications. As far as we know, there is no LB model for the two-component system of coupled sine-Gordon equations. The system has similar structure of the convection-diffusion system except for the second time derivative. We can define the first derivatives of macroscopic variable as the sum of distribution functions by the thought of the reference [44]. The main goal of this work is to extend the LB model to solve this two-component system of coupled sine-Gordon numerically by using the double mesoscopic Boltzmann equations. Through the C-E multiscale expansion, the governing nonlinear coupled evolution equations are recovered accurately from the double continuous Boltzmann equations. In order to compare the numerical solutions with the analytic ones, three test problems are taken into account. It is found that the numerical solutions are in accordance with the analytical ones. This demonstrates that the present model is an valid and flexible way for actual application.
The content of this paper is arranged as follows. Next section shows our LB model for the coupled sine-Gordon equations for the two-component system through the present model. Numerical validation is presented in Section 3. Finally, a brief summary is made.

Lattice Boltzmann Model
In the present model, the three-velocity lattice Bhatnagar-Gross-Krook (BGK) model is used. The directions of the particle discrete velocity are defined as e i , (i = 0, 1, 2): The LB equation with double distribution functions for u(x, t) and w(x, t) are given as follows (s = 1, 2): where f si (x, t) and f (0) si (x, t) refer to the distribution function and equilibrium distribution function, respectively. h si (x, t) is defined as an amending function, c is a constant to determine the viscous coefficient, ∆t is the time step, τ s the dimensionless single-relaxation-time which regulates the rate of approach to the equilibrium. The stability of the equation needs τ s > ∆t/2 [61].
Unlike the normal LB method, the first derivatives of macroscopic variables u(x, t) and w(x, t) are defined [44] as follows: Thus, the steady macroscopical quantities meet the following conservative conditions: Afterwards, through choosing appropriate local equilibrium distributions and amending functions, the corresponding macroscopic coupled system can be retrieved correctly.
Next, we will give the detailed derivation. Applying the Taylor expansion to the left-hand side of Equation (5) about the point x and t, we can obtain (s = 1, 2): By introducing the C-E multiscale expansions, we can expand the distribution function f si around si as follows: And f (k) si (x, t) (k = 1, 2, · · · ) are the non-equilibrium distribution functions, which satisfy the solvability conditions (s = 1, 2): Dividing both sides of Equation (8) by ∆t and substituting Equation (9) into Equation (8), we have: Comparing the two sides of Equation (11) and setting terms of order ε to each other, we have O(ε): Therefore: Comparing the two sides of Equation (11) and setting terms of order ε 2 to each other, we get O(ε 2 ): Substituting Equation (13) into Equation (14), we get: that is: Summing Equation (12) over i, we obtain: Summing Equation (16) over i, and using Equation (17), we obtain: According to the macroscopic equations, the local equilibrium distribution function f si (x, t) is required to satisfy the following relations: in terms of: Meanwhile, the source term h si satisfies: With Equation (19), Equation (17) becomes: and: With Equations (19), Equation (18) becomes: and: When (22) ×ε + (24) ×ε 2 is applied, the final equation is: When (23) ×ε +(25) ×ε 2 is applied, the final equation is: Meanwhile, from Equation (19), we can get the local equilibrium distribution functions si (x, t), (s = 1, 2, i = 0, 1, 2) as: From Equation (21), the amending functions h si (x, t), (s = 1, 2, i = 0, 1, 2) can be determined. For the sake of simplicity, only one case is presented here: and: In the simulation process, in order to get u(x, t) and w(x, t), we can apply the backward difference to the items ∂u(x, t) ∂t and ∂w(x, t) ∂t as: and: then using Equation (6), we get: and:

Numerical Simulation
In order to test the accuracy and efficiency of the present LB model, three initial and boundary value problems which have analytical solutions are simulated.
For the sake of numerical stability of the finite LB scheme, ∆x/∆t ≥ 1 is adopt in all simulations. Initially, the distribution functions f si (x, t) are set to equal f si (x, t). And the macroscopic variables u(x, t) and w(x, t) in Equation (1) are set to equal the initial conditions. The initial and boundary conditions of the test problems with analytical solutions are content with their analytical solutions. The non-equilibrium extrapolation scheme [62] is adopted to deal with the boundary condition.
Firstly, let us introduce symbols f n s,i,j = f si (x j , t n ), (s = 1, 2, i = 0, 1, 2), u n j = u(x j , t n ), w n j = w(x j , t n ), x j = j∆x, t n = n∆t, n is the nth layer time, j is the spatial grid. Then we can reformulate the LB Equation (5) by the classical finite difference notation: At time (n + 1)∆t, u and w are updated as follows: and: w n+1 The initial local equilibrium distribution functions f 0 s,i,j , (s = 1, 2, i = 0, 1, 2) are: The global relative error (GRE) is introduced to measure the present model's precision, and defined as follows: where u(x j , t), u * (x j , t) represent the numerical solution and analytical one, respectively. The summation is taken all grid points together. Next, numerical tests are performed for different initial conditions of the coupled sine-Gordon equations. It is found that the numerical solutions are in accordance with the analytical solutions over a relatively long period of time.

Example 1.
Consider the two-component system of coupled sine-Gordon equations in the region −5 ≤ x ≤ 5 given as: with the initial conditions: and the analytical solution for this problem is extracted from Ref. [56] by: where b = (1 − α 2 δ 2 )/(1 − δ 2 ). The boundary conditions conform to the analytical solution.
In the proceeding, we adopt α = 0.01, δ = 0.1, µ = 0.3, c 1 = c 2 = 1.0, τ 1 = τ 2 = ∆t. The computational region is fixed on I = [−5, 5]. The global relative errors (GRE) for the solutions u(x, t) and w(x, t) at t = 0.2 with different resolutions, from ∆x/∆t = 10 to 80, and the space grid N from 400 to 3200, are listed in Tables 1 and 2. From these two tables, we can see that GREs for u(x, t) are found to range from 9.9057 × 10 −7 to 6.3625 × 10 −5 , and GREs for w(x, t) are found to range from 1.8355 × 10 −2 to 1.8455 × 10 −2 . We can also see that when ∆x/∆t is larger, namely ∆t is relatively smaller, the global relative error of u(x, t) reduces with the first-order accuracy, while the global relative error of w(x, t) changes little. The accuracy of the macroscopic variable w(x, t) is not affected by the resolution in space and time. That is due to the treatment of the items ∂u(x, t)/∂t and ∂w(x, t)/∂t in Equations (31) and (32), has the first-order accuracy. According to Tables 1 and 2, we adopt N x × N t = 800 × 320 in consideration of both the computational accuracy and efficiency. It can be found that the numerical results in according with the analytical solutions, which are presented as the spatiotemporal evolution graph of the numerical (left) and analytical (right) solutions for comparison, see Figures 1 and 2. For clarity of contrast, we also present the two-dimensional visual comparisons of u(x, t) (left) and w(x, t) (right) at t = 0.2, see Figure 3. It is evident that the numerical results coincide with the analytical solutions.
In the proceeding, we adopt α = 0.01, δ = 0.05, µ = 0.2, c 1 = c 2 = 1.0, τ 1 = τ 2 = ∆t. The computational region is fixed within I = [−10, 10]. The GREs for the solutions u(x, t) and w(x, t) at t = 0.2 in different resolutions, from ∆x/∆t = 10 to 80 and the space grid N from 400 to 3200, are listed in Tables 3 and 4. From these two tables, we can see that the GREs for u(x, t) range from 1.0980 × 10 −6 to 7.0536 × 10 −5 , and the GREs for w(x, t) range from 2.8141 × 10 −2 to 3.2528 × 10 −2 . When ∆x/∆t is larger, namely ∆t is relatively small, the GREs of u(x, t) reduces with first-order accuracy, the GREs of u(x, t) change little. At the same time, we present the spatiotemporal evolution graph of the numerical (left) and analytical (right) solutions of u(x, t) and w(x, t) for comparison, see Figures 4 and 5. For clarity of contrast, we also present the two-dimensional visual comparisons of u(x, t) (left) and w(x, t) (right) at t = 0.2, see Figure 6. It can be found that the numerical results in according with the analytical solutions.

Example 3.
Consider the two-component system of coupled sine-Gordon equations in the region −10 ≤ x ≤ 10 given by: with the initial conditions: and the analytical solution for this problem is extracted from Reference [56] by: . The boundary conditions conform to the analytical solution.
In the proceeding, we adopt α = 1.6, δ = 2.0, b = 2.5, c 1 = c 2 = 1.0, τ 1 = τ 2 = ∆t. The computational region is fixed within I = [−10, 10]. The GREs for the solutions u(x, t) and w(x, t) at t = 0.2 in different resolutions, from ∆x/∆t = 10 to 80 and the space grid N from 400 to 3200, which are listed in Tables 5 and 6. From these two tables, we can see that the GREs for u(x, t) are found to range from 3.2169 × 10 −5 to 2.3832 × 10 −4 , and the GREs for w(x, t) are found to range from 2.9026 × 10 −3 to 2.9126 × 10 −3 . When ∆x/∆t is larger, namely ∆t is relatively small, the global relative error of u(x, t) reduces with first-order accuracy, the GREs of u(x, t) changes little. At the same time, we present the spatiotemporal evolution graph of the numerical (left) and analytical (right) solutions for comparison, see Figures 7 and 8. For clarity of contrast, we also present the two-dimensional visual comparisons of u(x, t) (left) and w(x, t) (right) at t = 0.2, see Figure 9. It can be found that the numerical results in according with the analytical solutions.   . Comparison between numerical and analytical solutions of u(x, t) (a) and w(x, t) (b) at t = 0.2 with N x × N t = 800 × 320 for Example 3. The blue circle symbol represents the numerical solution, and the solid red line represents the analytical solution given by Equation (48).

Conclusions
In conclusion, we have researched the application of the LB method for the solution of the two-component system of coupled sine-Gordon equations. By choosing the equilibrium distribution function and an amending function suitably, according to our proposed model, the governing evolution equations can be recovered accurately, in which the Chapman-Enskog multiscale expansion is employed. Numerical simulation for three test problems has been conducted to validate the LB model. The numerical results are in good agreement with the analytical solutions. While we take different initial conditions, we can get unique numerical solutions. We can also see that when ∆x/∆t is larger, ∆t is relatively smaller, and the global relative error of u(x, t) reduces with first-order accuracy; nevertheless, the global relative error of w(x, t) changes little. It is found that the accuracy of the macroscopic variable w(x, t) is not affected by the resolution in space and time. That is due to the treatment of the items ∂u(x, t)/∂t and ∂w(x, t)/∂t in Equations (31) and (32), which have first-order accuracy. For the purpose of attaining better computational accuracy and efficiency, the LB method for the test problems needs relatively small time step and space step. The present model can be developed to research more different types of the nonlinear system problems. There are still many problems worth studying to develop the present method, such as how to improve the accuracy and stability; we will continue these study in the near feature.

Conflicts of Interest:
The authors declare no conflict of interest.