A Fourth Order Energy Dissipative Scheme for a Trafﬁc Flow Model

: We propose, analyze and numerically validate a new energy dissipative scheme for the Ginzburg–Landau equation by using the invariant energy quadratization approach. First, the Ginzburg–Landau equation is transformed into an equivalent formulation which possesses the quadratic energy dissipation law. After the space-discretization of the Fourier pseudo-spectral method, the semi-discrete system is proved to be energy dissipative. Using diagonally implicit Runge–Kutta scheme, the semi-discrete system is integrated in the time direction. Then the presented full-discrete scheme preserves the energy dissipation, which is beneﬁcial to the numerical stability in long-time simulations. Several numerical experiments are provided to illustrate the effectiveness of the proposed scheme and verify the theoretical analysis.


Introduction
The Ginzburg-Landau equation (GLE) was first proposed by Ginzburg and Landau in 1950 and applied to model the electrodynamics, quantum mechanics and thermodynamic properties of superconductors. As a state transformation equation, the GLE is closely related to many other state transformation equations, such as Allen-Cahn equation and Chafee-Infante equation. As is well known, the GLE plays an important role in the study of state transformation and unstable wave theory. Moreover, it has been widely used in optical fiber communication. In this paper, we consider the one-dimensional GLE as a model for traffic flow [1] subject to the initial and boundary conditions u(x, 0) = u 0 (x), x ∈ Ω, u(a, t) = u(b, t), t ≥ 0, where ≥ 0 and Ω = [a, b] ⊆ R is the bounded domain. We denote the L 2 inner product by f , g = Ω f gdx, the L 2 norm f 2 = f , f where f (u) = u 4 . Then the GLE in [1] can be rewritten as Taking the inner product of the Equation (2) with m and the Equation (3) with u t , we obtain Obviously, the GLE satisfies the energy dissipation law.
In recent years, many researchers have conducted various works about the GLE. Celledoni [1] proved that the average vector field (AVF) method can keep the energy dissipation of the Equation (1). The time-dependent GLE for the car traffic model is derived by Nagatani [2] through the perturbation method. Gunzburger and Hou [3] solved a series of problems for a Ginzburg-Landau model of superconductivity. They proposed finite element approximations for the optimality system and derived corresponding error estimates. The local well-posedness of the initial-boundary value problem for the complex GLE in bounded domains is investigated by Kuroda and Takanori [4]. Meanwhile, there are many effective ways to solve the GLE, such as the discontinuous Galerkin method [5], the semi-implicit method [6] and the trial equation method [7]. In addition, there also exist relevant studies on the fractional GLE [8][9][10][11] and the stochastic GLE [12][13][14][15][16]. However, most of previous researches can only achieve second order accuracy in the time direction. Therefore, the main challenge issue is to construct a high order energy dissipative scheme for the GLE.
The invariant energy quadratization (IEQ) approach [17] is a new method developed in recent years by introducing a Lagrange multiplier and has been successfully applied to a variety of gradient flow models [18][19][20][21][22][23][24]. By introducing a Lagrange multiplier and converting the nonlinear part to a semi-implicit scheme, a new linear system consisting of two coupled equations is calculated at each time step. A time-discretization scheme for the linear systems is constructed by Yang [25], which is easy to implement and achieve unconditionally energy stability. Liu [26] applyed two novel and efficient modified techniques based on IEQ approach to deal with nonlinear terms in gradient flows which succeeded in finding suitable positive preserving functions. Two linear and second order schemes that combine the IEQ approach with the stabilization technique are developed by Xu [27], where several extra linear stabilization terms are added and they are crucial to suppress the non-physical spatial oscillations caused by the strong anisotropy. Gong [28] presented a novel class of arbitrarily high-order and unconditionally energy-stable algorithms for gradient flow models by combining the IEQ approach and a specific class of Runge-Kutta (RK) methods. Since the solution to the Ginzburg-Landau equation shows instability very soon and this high order implicit scheme can unconditionally keep Ginzburg-Landau equation energy dissipation, we adopt this method to solve the Ginzburg-Landau equation. Ref. [28] combines IEQ approach and RK method and apply them to Cahn-Hilliard equation and Allen-Cahn equation to get the arbitrarily high-order schemes, while we apply it to the Ginzburg-Landau equation to construct a fourth order scheme in this paper.
The rest of this paper is organized as follows. In Section 2, the IEQ approach is introduced and applied to the GLE. Then the GLE is reformulated as a new system which admits a quadratic energy dissipation law. In Section 3, the Fourier pseudo-spectral method is used to obtain the semi-discrete the system of the GLE. After that, the implicit fourth order Runge-Kutta (RK4) method is applied for the time discretization to structure a full-discrete scheme. It is proved that the scheme we proposed is energy dissipative. In Section 4, the stability of the proposed scheme is verified by numerical experiments. We make a summary of this work in Section 5.

Invariant Energy Quadratization Approach
In this section, we present a new scheme for the GLE to preserve the energy dissipation law. By using the IEQ approach, we transform the Hamlitonian functional into a quadratic functional. We first introduce a Lagrange multiplier r = f (u), then the GLE is transformed into Theorem 1. The system above preserves the energy dissipation law

Proof.
The system (4) can be rewritten as By taking the inner products of the Equations (5)- (7) with m, u t and r 2 , respectively, we have the following equalities By substituting Equations (9) and (10) into Equation (8), we obtain

Spatial Discretization by Using Fourier Pseudo-Spectral Method
In this part, we choose the Fourier pseudo-spectral method to discretize the system (4). It is proved that the semi-discrete scheme for the GLE preserves the semi-discrete energy dissipation law.
We primarily divide the domain Ω = [a, b] into N equal intervals, where N is an even number. Afterwards, denote h = b−a N as the spatial step and x j := a + jh(j = 0, 1, · · · , N − 1) as the grid points. The kth-order partial differential optrator ∂ k x can be approximated by the Fourier spectral differential matrix D k where µ = 2π b−a , j, l = 0, 1, · · · , N − 1. Clearly, D 1 is an N × N skew-symmetric matrix. By applying Fourier pseudo-spectral method to the system (4), we acquire the semi-discrete scheme Theorem 2. The system (11) preserves the semi-discrete energy dissipation law Proof. The system (11) can be rewritten as follows Taking the inner products of the system (12) with m, u t and r 2 , respectively, we obtain the following semi-discrete energy dissipation law It is known that the Equation (13) is equivalent to where is a symmetric and constant matrix.

Full Discretization by Using Implicit High Order Runge-Kutta Method
In previous studies, many methods have been proposed to solve energy dissipation system, such as the second order backward differentiation formula (BDF2), the AVF method and the backward Euler method. However, these methods can only structure first or second order energy dissipative schemes in the time direction meanwhile these schemes always require a lot of iterative computations, which not only affect the efficiency, but also bring cumulative errors in the calculation processes. Therefore, we propose a new scheme to acquire high order energy dissipative schemes for the GLE.
Let b i , a ij , (i, j = 1, 2, · · · , s) be real numbers and c i = ∑ s j=1 a ij . For given z n = z(t n ), when an s-stage RK method is applied to the Equation (1), the following intermediate values are first calculated by where F j = F(z j ). Then the high order energy dissipation scheme of the GLE is updated via where A ∈ R s,s , b ∈ R s and c = AI with I = (1, 1, · · · , 1) T ∈ R s .
In consequence, we choose RK(4,4) as Definition 1. Let S be a symmetric, constant matrix. For a dissipative system, the RK method defined by the Butcher array is said to be S-dissipative if applid to the new system (15) Sz n+1 , z n+1 ≤ Sz n , z n for all n ≥ 0 and for all ∆t ≥ 0.
are nonnegative.
Proof. Above all, replace z n+1 by the new system (15), that can be derived as Substituting z n = z i − ∆t ∑ s j=1 a ij F j into the first term on the right hand side of the Equation (17), we get Whereafter, it can be arranged as Noticing that Sz i , F j ≤ 0 and m ij ≥ 0, we have Sz n+1 , z n+1 ≤ Sz n , z n . This completes the proof.
In this paper, we only consider four stage and fourth order diagonally implicit RK scheme, which belongs to S-dissipative RK scheme can be extended to higher order RK method without any necessary difficulties.

Numerical Simulation
In this section, we carry out two numerical experiments by applying the implicit high order energy dissipative scheme. The stability and the convergence of the new scheme (15) are verified by these numerical experiments. Meanwhile, we compare the proposed scheme (15) with the AVF method in energy errors and accuracy. The numerical iteration formula of the AVF method is where E = D 1 + D 2 . In the following text, the Newton iteration is mainly used to solve the full-discrete system.

Example 1.
We utilize the initial condition of single solitary wave given by to observe how the coefficient affects the energy dissipation. First, we set M 1 = 0.3, the number of grid points N = 200 and the time step ∆t = 0.001. Letting ε = 0, 0.00025, 0.0005, 0.00075, 0.001 respectively, the energy variations and energy errors of the proposed scheme (15) are presented in Figure 1. Based on Figure 1a, we find the proposed scheme (15) can preserve the discrete energy conservation exactly when = 0, and the energy dissipations of the GLE can be observed when > 0. Obviously, with bigger , the energy dissipates faster. Taking the discrete energy obtained from the scheme (15) with ∆t = 0.0001 and N = 200 as the reference discrete energy, we plot the energy errors (∆t = 0.001) with various ε in Figure 1b. It shows that the proposed new scheme (15) can preserve the energy dissipation of the GLE. Meanwhile, Figure 1b indicates that the energy dissipation is a little bit influenced by the time step size. The waveforms, energy variations and energy errors obtained from the scheme (15) and the scheme (18) with ε = 0.001, N = 100, ∆t = 0.001 are pictured in Figure 2. The waveforms obtained from these two schemes at t = 0.1 are almost the same (Figure 2a). Figure 2b shows that they both can keep the energy reducing steadily. According to Figure 2c, comparatively speaking, the energy errors obtained by the proposed scheme (15) is even smaller than that by the scheme (18 Table 1, the convergence order of the proposed scheme (15) in time is 4 and the scheme (18) only converges with second-order accuracy. For the space accuracy test, the scheme (15) and the scheme (18) both apply the Fourier pseudo-spectral method in space-discretization, thus the space convergence orders of them decrease in an exponential rate in Table 2 (the numerical solution with N = 560 and ∆t = 0.001 is chosen as the reference solution).
By setting M 1 = M 2 = −1 and ε = 0.001, we plot the waveform of double solitary wave with N = 200 and ∆t = 0.001 in Figure 3a. The energy dissipation property of the GLE with various time steps can be verified in Figure 3b. It shows the new scheme (15) can keep energy dissipation at a large time step. Taking the discrete energy obtained from the new scheme (15) with ∆t = 0.0001 as the reference discrete energy, the energy error in Figure 3c is only about O(10 −7 ). Then using the same grid, we consider the situation of M 1 = −0.2 and M 2 = −2 and plot the energy variation and the energy error in Figure 4. As we can see, the proposed scheme (15) still can maintain the energy dissipation property with various time steps.    Energy error

Conclusions
In this paper, we have developed a stable and high order energy dissipative scheme for solving the GLE. By using the IEQ approach, the GLE is downgraded from fourth order to second order. The full-discrete system of GLE is obtained by using the Fourier pseudo-spectral method in the spatial direction and the diagonally implicit Runge-Kutta method in the time direction. The new scheme (15) we proposed has three main advantages. First of all, it is well-known that the implicit methods enjoy better stability in the simulation of stiff equations. Besides, it has high calculation efficiency due to the application of the diagonally implicit Runge-Kutta method. Last but not least, the high order scheme in time direction can be constructed and the energy dissipation can be preserved. The proposed implicit fourth order scheme can maintain the energy dissipation of the GLE theoretically. Moreover, the numerical results of single solitary wave and double solitary wave have verified the stability and accuracy of the proposed scheme. Compared with the scheme (18), the proposed scheme converges with fourth order in the time direction and can be generalized to construct higher order energy dissipative scheme naturally.

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