A Parallel Spectral Element Method for Fractional Lorenz System

We provide a parallel spectral element method for the fractional Lorenz system numerically. The detailed construction and implementation of the method are presented. Thanks to the spectral accuracy of the presented method, the storage requirement due to the “global time dependence” can be considerably relaxed. Also, the parareal method combining spectral element method reduces the computing time greatly. Finally, we have tested the chaotic behaviors of fractional Lorenz system. Our numerical results are in excellent agreement with the results from other numerical methods.


Introduction
In 1963, Lorenz [1] attempted to set up a (much simpler and easier-to-analyze) system of differential equations that would explain some of the unpredictable behavior of the weather.The classical Lorenz system is where all three parameters, , , , are assumed to be positive, and moreover,  > +1.Lorenz made an important discovery that all nonequilibrium solutions tend eventually to the same complicated set, the Lorenz attractor.The behavior of the chaotic system is one of the most important subjects in modern mathematics.
In [2], I. Grigorenko and E. Grigorenko generalized the Lorenz system to the fractional case by replacing the usual derivatives by the fractional one and found some interesting chaotic behaviors of fractional Lorenz system.This generation is natural.As we know, fractional derivatives play a more and more important role in modern sciences.They appear in the investigation of transport dynamics in complex systems which are characterized by the anomalous diffusion and nonexponential patterns [3].It has been shown that many important physical and biological systems can be described by fractional differential equations (see [3][4][5][6], etc.).This has led to some investigations on the theoretical analysis and scientifical computation of the fractional differential equations in recent years (for details see [7][8][9][10][11][12][13] and the references therein).However, the fact that the fractional derivative is a nonlocal operator makes the results of fractional calculus in many respects far from complete.
First, let us recall that, for 0 <  < 1 and () defined on the interval [, ], the (Caputo) fractional derivative is defined as where, Γ(⋅) is the Gamma function.And for simplicity, we denote () = 1/Γ(1 − ).Thus, the fractional generalization of the Lorenz system is where 0 < , ,  ≤ 1, and  ≥ 1 is an integer.In [2], the authors pointed out that the system can have an effective noninteger dimension Σ defined as the sum of the orders of all involved derivatives.They found that the system with Σ < 3 can also exhibit chaotic behavior.Furthermore, a striking finding is that there is a critical value of the effective dimension Σ cr , under which the system undergoes a transition from chaotic daynamics to regular one.The numerical scheme of [2] is based on the linearization of system (3), on each step of integration and iterative application of the exact solution of linear fractional differential equations.Also, homotopy methods to solve the (fractional) Lorenz system were proposed in [14,15].It has been known that the feature of the fractional derivatives makes the design of accurate and fast methods difficult.Unlike the integer derivatives, which are local in the sense that the derivative of a function at a certain point depends only on the function in the vicinity of this point, presence of the integral in the noninteger order derivatives makes the problem global.This means that the solution at a time   depends on the solutions at all previous time levels  <   .The fact that all previous solutions have to be saved to compute the solution at the current time level would make the storage very expensive if low-order methods are employed for temporal discretization.Thus, it is very desirable to use high-order methods for efficient computations of the numerical solution of this kind of problem.Among the high-order methods, spectral methods are known to be very useful for their exponential rate of convergence, which is also demonstrated for solving fractional differential equations [8,9].Thanks to the spectral accuracy, the storage requirement due to the "global time dependence" can be considerably relaxed.However, a drawback of the spectral method is also wellknown; that is, the matrix associated with the spectral method is full and the computational cost grows more quickly than that of a low-order method.Thus for long integration it is desirable to combine the spectral method with domain decomposition techniques to avoid using a single high-degree polynomial (see [16] for details).This motivates us to consider the parallel spectral element method, since, as compared to a low-order method, this method needs fewer grid points to produce highly accurate solutions, and, furthermore, the parallel method reduces the computing time greatly.
This paper is organized as follows.In Section 2, we construct and implement the parallel spectral element schema for the numerical solution of the fractional Lorenz system.In Section 3, we have tested the chaotic behavior of fractional Lorenz system using our numerical method.The numerical results are in excellent agreement with the results from other numerical methods.Finally, Section 4 includes some concluding remarks.
The main idea of the parallel method consists of two steps.First, design two families of discrete solvers, fine solver F  and coarse solver G  , to approximate   .The fine approximate solver F  is defined by the spectral collocation method with polynomial degree of : with the polynomial U  = (  ,   ,   ) of degree .The coarse solver operator G  is defined in the same way, using polynomials of a lower degree Ñ < .Second, set up the iteration: with the initial value where subscript  refers to the element number, subscript  refers to the iteration number, and the polynomial U   = (   ,    ,    ) is the approximate solution of X  at the time interval   .The parallel scheme can be seen as the predictorcorrector method, with the predictor G  (U  1 , . . ., U  −1 ) and the corrector The key point of the parallel method is that the corrector C  can be computed in parallel for 1 ≤  ≤ .The predictor term G  (U  1 , . . ., U  −1 ) in ( 7) must be computed sequentially, but this is relatively cheap if Ñ is small compared to .Thus, provided the iteration converges rapidly, we can use multiple processors to obtain the high-accuracy solution in a small multiple of the time needed to compute a lowaccuracy solution.
Now we begin to describe the fine approximate solver F  .We first define P  as the polynomial space of degree less than or equal to .We denote by  ,  (),  ∈ [−1, 1], the Jacobi polynomial of degree  with respect to the weight function  , () = (1 − )  (1 + )  , −1 < ,  < 1.Let  ,  be the points of the Gauss-Labatto-Jacobi (GJ) quadrature formula, defined by arranged by increasing order: .The associated weights of the GLJ quadrature formula are denoted by  ,  , 0 ≤  ≤ .Particularly, let    be the Gauss-Labatto-Legendre (GLL) points on the interval   .

Now we turn to compute ∫
Thus, applying (18), (20), and (22) yields Similarly, we can compute the matrices  2 , and  3 , .Equation ( 13) is a nonlinear nonsymmetric system.We use the Jacobian-free Newton-Krylov methods to solve it.Let U  = (u  , k  , w  )  ; we rewrite (13) in the following form: The iteration can be described as follows: (i) The initial guess U (0)  is given.(ii) Using BI-Conjugate-Gradient-stab (BICGstab) method solve the following linear system: where U = (u, v, w)  ,  is the iteration index, and J   (U () ) is the Jacobian matrix of function   at U ()  .We need not compute the matrix J   (U ()  ).In fact, only the matrix-vector product is needed in the BICGstab method.Let   (U ()   , U) be the Fréchet derivative of the function   at U ()   with respect to the direction U.Consider (iii) Update U ()  by where  is the correction factor, which is selected such that       (U (+1) where  = 1/10 is a scale parameter; we test  = 1, 1/4, 1/16, . .., until (28) is satisfied.
(iv) Stop the iteration, if
The phase portraits of , ,  are shown in Figures 1 and 2. From these figures, it may be concluded that, for  =  =  = 0.98 up to  =  =  = 1, the system is chaotic, while for  =  =  < 0.97, the system turns to be regular and the solution converges to one of the two equilibrium points  ± .This also coincides with other numerical results.
Secondly, we take  = 2, while other parameters ,  and  are given above.Furthermore, we take  = 3, while other parameters ,  and  are given above.The equilibrium points of (3) are the origin, and The results are similar to the case of  = 1.However, the critical efficient dimension Σ cr is much smaller than the case  = 1.From Figure 3, it may be concluded that, for  =  =  = 0.94 up to  =  =  = 1, the system is chaotic, while for  =  =  < 0.93, the system turns to be regular and the solution converges to one of the two equilibrium points  ± .Furthermore, if we take  = 3, while other parameter as given above.The equilibrium points are the origin, and We can see from Figure 4 that the critical efficient dimension Σ cr is much smaller than the case  = 2.

Comparison with Other Methods.
In Table 1, we compare the spectral method with classical finite difference method and finite element method with respect to convergence rate, computational complexity, and storage.Spectral method is well-known for its exponential convergency, that is,  − , where  is the freedom of the discretized space.The convergence rate is slower than  −2 for finite element method [7,12] and slower than  − ,  ≤ 4 for finite difference method [10,[17][18][19][20].
For partial differential equations of integral order, the final matrix is full for spectral method while sparse for both finite difference method and finite element method.Compared with spectral method, sparse matrix is the main advantage of finite difference method and finite element method for solving integral differential equations.However, for the fractional differential equations, the final matrix for all three methods is full.The presence of the integral in the noninteger order derivatives means that the solution at a time   depends on the solutions at all previous time levels  <   .This fact makes the final matrix full.So the computation complexity is about ( 3 ) with ( 2 ) for matrix-vector multiple and () for iteration steps.
Furthermore, all previous solutions have to be saved to compute the solution at the current time level which would make the storage very expensive if low-order methods are employed for temporal discretization.For example, we consider fractional ordinary differential equations: ), and (10 5 ) mounts of storage are used to approximate an exact solution by spectral method, finite difference method [10] ( = 2), and finite element method [12] ( = 2), respectively.In view of these three advantages, spectral method is the national choice for the fractional differential equation.

Efficiency of the Parallel Method.
The parallelism efficiency is demonstrated by a cost comparison between the parallel in time scheme and the classical sequential scheme based on the corresponding fine mesh.
First, the classical sequential scheme based on the fine mesh consists of solving the problems F  (,  1 ,  2 , . . .,  −1 ) consecutively for 1 ≤  ≤ .The computational complexity is equal to the sum of all the elementary operations: If we implement the scheme in a parallel architecture with enough processors, then the total computational time corresponds to the cost to solve a sequential set of  coarse subproblems and a fine subproblem in a single processor.If  is the number of iterations required to achieve the desired convergence of the parallel in time algorithm, then the total computational complexity is Comparing (34) with (35), we obtain a speed up (i.e., the percentage with respect to the sequential scheme) close to

Table 1 :
Comparison with FDM and FEM.