A Parallel Spectral Element Method for Fractional Lorenz System

Abstract

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.

1. 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, . 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–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–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 and defined on the interval , the (Caputo) fractional derivative is defined as where, is the Gamma function. And for simplicity, we denote . Thus, the fractional generalization of the Lorenz system is where , and 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 can also exhibit chaotic behavior. Furthermore, a striking finding is that there is a critical value of the effective dimension , 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 well-known; 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.

2. The Implementation of the Parallel Spectral Element Method

In this section, we implement the parallel spectral element method for fractional Lorenz system (3). First, we divide the interval into subintervals , . The length of each subinterval is , and , . According to this partition, (3) is divided into the following system of equations: with , , and , . Let ; then the solution of (4) can be expressed by the solvers as follows: In most cases is not realizable and needs to be approximated.

The main idea of the parallel method consists of two steps. First, design two families of discrete solvers, fine solver and coarse solver , to approximate . The fine approximate solver is defined by the spectral collocation method with polynomial degree of : with the polynomial of degree . The coarse solver operator 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 is the approximate solution of at the time interval . The parallel scheme can be seen as the predictor-corrector method, with the predictor and the corrector . The key point of the parallel method is that the corrector can be computed in parallel for . The predictor term 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 low-accuracy solution.

Now we begin to describe the fine approximate solver . We first define as the polynomial space of degree less than or equal to . We denote by , , the Jacobi polynomial of degree with respect to the weight function , . 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 , . Particularly, let be the Gauss-Labatto-Legendre (GLL) points on the interval .

Thus, the spectral collocation method for fractional Lorenz system (3) on each subinterval is to find , such that for all

For the nonlinear problem (10), we choose to use the Lagrangian polynomials as a basis of the approximation spaces. Let be the Lagrangian polynomials associated with the points . That is, such that , where denotes the Kronecker function. It is seen that the set forms a basis of : By expressing , , and and all nonlinear term in this basis, We arrive at the matrix statement of (10): where , , and are the unknown vectors. is acting on two vectors . The matrices , , , , , and are matrices. They are computed as follows:

We are now led to compute the integrals and one by one.

For the first one, we employ Gauss-type quadrature to calculate it. It is well-known that the numerical quadrature is exact for all functions such that . As a result, by the linear variable transformation it holds since for all .

For the second one, it is important to point out that the Gauss-type quadrature fails if approach . Here we design a new technique to calculate this integral precisely. We first rewrite the integral into the form by the linear transformation where is the Lagrangian polynomials based on the GLL points , , and . Next, we express by the Legendre basis: where .

Now we turn to compute . By employing the property and integrating by parts, we get Thus, applying (18), (20), and (22) yields Similarly, we can compute the matrices and .

Equation (13) is a nonlinear nonsymmetric system. We use the Jacobian-free Newton-Krylov methods to solve it. Let ; we rewrite (13) in the following form: for ,

The iteration can be described as follows:(i)The initial guess is given.(ii)Using BI-Conjugate-Gradient-stab (BICGstab) method solve the following linear system:  where , is the iteration index, and is the Jacobian matrix of function at . We need not compute the matrix . In fact, only the matrix-vector product is needed in the BICGstab method. Let be the Fréchet derivative of the function at with respect to the direction . Consider (iii)Update by  where is the correction factor, which is selected such that  where is a scale parameter; we test , until (28) is satisfied.(iv)Stop the iteration, if

3. Numerical Results

For a comparison with [2, 14], we take the initial values , , and . In our numerical experiments, we take , , , and .

First, we consider the case , , , and . Note that, in this case, the equilibrium points of (3) are the origin, and As pointed in [2], the solution can exhibit chaotic behavior if the effective dimension is greater than the critical dimension . In our numerical experiments, we take and and . The phase portraits of are shown in Figures 1 and 2. From these figures, it may be concluded that, for up to , the system is chaotic, while for , the system turns to be regular and the solution converges to one of the two equilibrium points . This also coincides with other numerical results.

(a)

(b)

Figure 1 

The phase portraits of -- for fractional Lorenz system (left) and (right) for the case .

Figure 2 

The phase portraits of -- for fractional Lorenz system for the case .

Secondly, we take , while other parameters , and are given above. Furthermore, we take , while other parameters , and are given above. The equilibrium points of (3) are the origin, and The results are similar to the case of . However, the critical efficient dimension is much smaller than the case . From Figure 3, it may be concluded that, for up to , the system is chaotic, while for , the system turns to be regular and the solution converges to one of the two equilibrium points .

(a)

(b)

Figure 3 

The phase portraits of -- for fractional Lorenz systems (left) and (right) for the case .

Furthermore, if we take , while other parameter as given above. The equilibrium points are the origin, and We can see from Figure 4 that the critical efficient dimension is much smaller than the case .

(a)

(b)

Figure 4 

The phase portraits of -- for fractional Lorenz systems (left) and (right) for the case .

4. Efficiency of the Method

4.1. 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.

Table 1 

Comparison with FDM and FEM.

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 for finite element method [7, 12] and slower than for finite difference method [10, 17–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 with 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: with , . If precision is desired, , then , , and mounts of storage are used to approximate an exact solution by spectral method, finite difference method [10] (), and finite element method [12] (), respectively. In view of these three advantages, spectral method is the national choice for the fractional differential equation.

4.2. 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 consecutively for . 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

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The author is very thankful to the anonymous referees, who carefully read through the paper and made helpful suggestions that have significantly improved the presentation of this paper. This work was partially supported by National NSF of China (11226081) and NSF of Fujian Province (2013J05003).