About this Journal Submit a Manuscript Table of Contents
International Journal of Differential Equations
Volume 2012 (2012), Article ID 495202, 19 pages
http://dx.doi.org/10.1155/2012/495202
Research Article

A Time-Space Collocation Spectral Approximation for a Class of Time Fractional Differential Equations

Department of Mathematics, School of Sciences, South China University of Technology, Guangzhou 510641, China

Received 22 May 2012; Revised 26 July 2012; Accepted 29 July 2012

Academic Editor: Fawang Liu

Copyright © 2012 Fenghui Huang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A numerical scheme is presented for a class of time fractional differential equations with Dirichlet's and Neumann's boundary conditions. The model solution is discretized in time and space with a spectral expansion of Lagrange interpolation polynomial. Numerical results demonstrate the spectral accuracy and efficiency of the collocation spectral method. The technique not only is easy to implement but also can be easily applied to multidimensional problems.

1. Introduction

Fractional differential equations have attracted in recent years considerable interest because of their ability to model complex phenomena. For example, fractional derivatives have been used successfully to model frequency-dependent damping behavior of many viscoelastic materials. They are also used in modeling of many chemical processes, mathematical biology, and many other problems in engineering. Related equations of importance are fractional diffusion equations, the fractional advection-diffusion equation for anomalous diffusion with sources and sinks, and the fractional Fokker-Planck equation for anomalous diffusion in an external field, and so forth.

In this paper, we consider the following time fractional differential equation (TFDE) 𝐷𝜇𝑡𝑢(𝑥,𝑡)=𝜆2𝑢(𝑥,𝑡)𝜈𝜕𝑢(𝑥,𝑡)𝜕𝜕𝑥+𝐷2𝑢(𝑥,𝑡)𝜕𝑥2+𝑓(𝑥,𝑡)=𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡),𝑎<𝑥<𝑏,𝑡>0,(1.1) where 𝑢(𝑥,𝑡)=𝜆2𝑢(𝑥,𝑡)(𝜈(𝜕𝑢(𝑥,𝑡)/𝜕𝑥))+(𝐷(𝜕2𝑢(𝑥,𝑡)/𝜕𝑥2)) is a linear differential operator. 𝜆,𝜈0, 𝐷>0 are given constants, 0<𝜇1, 𝑓(𝑥,𝑡) is a given continuous function, 𝐷𝜇𝑡𝑢(𝑥,𝑡) is a time fractional derivative which is defined in the Caputo sense 𝐷𝛼𝑡𝜕𝑓(𝑡)=𝑚𝑓(𝑡)𝜕𝑡𝑚1,𝛼=𝑚,Γ(𝑚𝛼)𝑡0(𝑡𝜏)𝑚𝛼1𝜕𝑚𝑓(𝜏)𝑑𝜏𝑚𝑑𝜏,𝑚1<𝛼<𝑚.(1.2) The use of Caputo derivative in the above equation is partly because of the convenience to specify the initial conditions [1].

The TFDE (1.1) includes a few special cases: time fractional diffusion equation, time fractional reaction-diffusion equation, time fractional advection-diffusion equations, and their respective corresponding integer-order partial differential equations.

There are many analytical techniques for dealing with the TFDE, such as integral transformation method (including Laplace’s transform, Fourier’s transform, and Mellin’s transform) [15], operational calculus method [6], Adomian decomposition method [7], iteration method and series method [8], and the method of separating variables [9].

One of the key issues with numerical solution of the TFDE (1.1) is design of efficient numerical schemes for time fractional derivative. Until now, most numerical algorithms have relied on the finite difference (FD) methods to discretize the fractional derivatives, and the numerical accuracy always dependent on the order of the fractional derivatives. On the other hand, those FD methods have been generally limited to simple cases (low dimension or small integration) and are very difficult to improve the numerical accuracy [1014]. Some numerical schemes using low-order finite elements (FE) have also been proposed [1517]. The fractional derivatives are defined using integrals, so they are nonlocal operators. This nonlocal property means that the next state of a system not only depends on its current state but also on its historical states starting from the initial time. This nonlocal property is good for modeling reality, but they require a large number of operations and a large memory storage capacity when discretized with low-order FD and FE schemes. From this point, the “global method”—the nonlocal methods, like the spectral method—is well suited to discretize the nonlocal operators like fractional-order derivatives. These methods naturally take the global behavior of the solution into account and thus do not result in an extra computational cost when moving from an integer order to a fractional-order model. For example, Hanert has proposed a pseudospectral method based on Chebyshev basis functions in space and Mittag-Leffler basis functions in time to discretize the time-space fractional diffusion equation [18, 19]. Li and Xu have proposed a Galerkin spectral method based on Lagrangian basis functions in space and Jacobi basis functions in time for time fractional diffusion equation [20].

In this paper, we propose a time-space collocation spectral method to discretize the TFDEs (1.1), which is easier to implement and apply to multidimensional problems than the existing Galerkin spectral. Another advantage of the present scheme is that the method can easily handle all kinds of boundary conditions.

2. Analytical Solution of the TFDE in a Bounded Domain

In this section, we present some analytical solutions of the TFDE which will be found helpful in the comprehension of the nature of such a problem.

We consider the TFED (1.1) with initial condition 𝑢(𝑥,0)=𝜙(𝑥),𝑥(𝑎,𝑏),(2.1) and Dirichlet boundary conditions 𝑢(𝑎,𝑡)=𝜑1(𝑡),𝑢(𝑏,𝑡)=𝜑2(𝑡),𝑡(0,𝑇),(2.2) or Neumann boundary conditions 𝑢𝑥(𝑎,𝑡)=𝜑1(𝑡),𝑢𝑥(𝑏,𝑡)=𝜑2(𝑡),𝑡(0,𝑇).(2.3)

For the case that 𝑓0 and 𝑎=0,𝜑1(𝑡)=𝜑2(𝑡)=0, by applying the finite sine (cosine) and Laplace transforms to (1.1) with initial condition (2.1), the analytical solutions for the problem can be obtained [5] as 𝑢(𝑥,𝑡)=2𝑒𝜈𝑥/2𝐷𝑏𝐷𝑛=1𝐸𝜇𝜈24𝐷+𝜆2+𝑛𝜋𝑏2𝑡𝜇sin𝑛𝜋𝑥𝑏𝐷×𝑏𝐷0𝜙(𝑦)sin𝑛𝜋𝑦𝑏𝐷𝑒𝜈𝑦/2𝐷𝑑𝑦,(2.4) for homogeneous Dirichlet boundary conditions, and 𝑢(𝑥,𝑡)=2𝑒𝜈𝑥/2𝐷𝑏𝐷𝑛=1𝐸𝜇𝜈24𝐷+𝜆2+𝑛𝜋𝑏2𝑡𝜇cos𝑛𝜋𝑥𝑏𝐷×𝑏𝐷0𝜙(𝑦)cos𝑛𝜋𝑦𝑏𝐷𝑒𝜈𝑦/2𝐷𝑑𝑦,(2.5) for homogeneous Neumann boundary conditions. Where 𝐸𝛼(𝑧) denotes a one-parameter Mittag-Leffler function which is defined by the series expansion 𝐸𝛼(𝑧)=𝑘=0𝑧𝑘Γ(𝛼𝑘+1),𝑧,(𝛼>0).(2.6)

Obviously, if we fix the variable 𝑥=𝑥, that is, 𝑢(,𝑡) is a function of the variable 𝑡, we can see the solution 𝑢(,𝑡) is not smooth on [0,𝑇]. According to (2.4) and (2.5), its first derivative behaves like 𝑢(,𝑡)𝑡𝜇1 and the high-order derivative behaves like 𝑢(𝑚)(,𝑡)𝑡𝜇𝑚 near 𝑡=0+.

3. Collocation Spectral Method

First, we give the properties of the Caputo fractional derivative [1] as 𝐽𝛽𝐷𝛽𝑔(𝑡)=𝑔(𝑡)𝑛1𝑘=0𝑡𝑔(0)𝑘𝑘!,0𝑛1<𝛽<𝑛,(3.1) where 𝐽𝛽 is the Riemann-Liouville fractional integral of order 𝛽 which is defined by 𝐽𝛽1𝑔(𝑡)=Γ(𝛽)𝑡0(𝑡𝜏)𝛽1𝑔(𝜏)𝑑𝜏.(3.2)

By the above properties, we can transform the initial value problem (1.1) into the following Volterra integral equation equivalently: 1𝑢(𝑥,𝑡)𝑢(𝑥,0)=Γ(𝜇)𝑡0𝑢(𝑥,𝜏)(𝑡𝜏)1𝜇1𝑑𝜏+Γ(𝜇)𝑡0𝑓(𝑥,𝜏)(𝑡𝜏)1𝜇𝑑𝜏.(3.3)

For the singular behavior of the exact solution near 𝑡=0+ which we have mentioned in the special case (the exact solution (2.4) or (2.5) behaves 𝑢(,𝑡)𝑡𝜇1 near 𝑡=0+), the direct application of the spectral methods is difficult. To overcome this difficulty, we use the technique in [21], that is, applying the transformation 𝑢(𝑥,𝑡)=𝑡1𝜇[],𝑢(𝑥,𝑡)𝑢(𝑥,0)(3.4) to make the solution smooth. Then (3.3) is transformed to the equation 𝑢(𝑥,𝑡)=𝑡𝑓(𝑥,𝑡)+1𝜇Γ(𝜇)𝑡0𝑢(𝑥,𝜏)(𝑡𝜏)1𝜇𝜏1𝜇𝑑𝜏,(3.5) where 𝑡𝑓(𝑥,𝑡)=1𝜇Γ(𝜇)𝑡0𝑓(𝑥,𝜏)𝑑𝜏(𝑡𝜏)1𝜇+𝑢(𝑥,0)Γ(1+𝜇)𝑡.(3.6)

To apply the theory of orthogonal polynomials, we set 𝑇𝑡=2[]𝑇(1+𝑦),𝑦1,1;𝜏=2[],(1+𝑠),𝑠1,𝑦(3.7) then the singular problems (3.5) can be rewritten as 𝑣(𝑥,𝑦)=𝑔(𝑥,𝑦)+(𝑇(1+𝑦)/2)1𝜇𝑇Γ(𝜇)22𝜇1𝑦1𝑣(𝑥,𝑠)(𝑦𝑠)1𝜇(1+𝑠)1𝜇𝑑𝑠,(3.8) where 𝑦[1,1], and 𝑣(𝑥,𝑦)=𝑢𝑇𝑥,2(1+𝑦),𝑔(𝑥,𝑦)=𝑓𝑇𝑥,2.(1+𝑦)(3.9)

For the collocation methods, (3.8) holds at the Gauss-Lobatto collocation points {𝑥𝑖}𝑁𝑖=0 and Jacobi collocation points {𝑦𝑗}𝑁𝑗=0 with Jacobi weight functions 𝜔(𝑦)=(1𝑦2)𝜇1 on [1,1], namely, 𝑣𝑥𝑖,𝑦𝑗𝑥=𝑔𝑖,𝑦𝑗+𝑇1+𝑦𝑗/21𝜇𝑇Γ(𝜇)22𝜇1𝑦𝑗1𝑥𝑣𝑖,𝑠𝑑𝑠𝑦𝑗𝑠1𝜇(1+𝑠)1𝜇.(3.10)

By using the following variable change: 𝑠=𝑠𝑗(𝜃)=1+𝑦𝑗2𝑦𝜃+𝑗12[],,𝜃1,1(3.11) we can rewrite (3.10) as follows: 𝑣𝑥𝑖,𝑦𝑗𝑥=𝑔𝑖,𝑦𝑗+𝑇1+𝑦𝑗/21𝜇Γ𝑇(𝜇)1+𝑦𝑗/222𝜇111𝑥𝑣𝑖,𝑠𝑗(𝜃)1𝜃21𝜇𝑥𝑑𝜃=𝑔𝑖,𝑦𝑗+𝑇1+𝑦𝑗/2𝜇22𝜇1𝑥Γ(𝜇)𝑣𝑖,𝑠𝑗(),1𝜔.(3.12)

We first use 𝑣𝑗𝑖, 0𝑖𝑁; 0𝑗𝑀 to indicate the approximate value for 𝑣(𝑥𝑖,𝑦𝑗), then we can use 𝑣𝑀𝑁(𝑥,𝑦)=𝑁𝑀𝑛=0𝑚=0𝑣𝑚𝑛𝑛(𝑥)𝐹𝑚(𝑦),(3.13) to approximate the function 𝑣(𝑥,𝑦), where 𝑛(𝑥) is the 𝑛th Lagrange interpolation polynomial associated with the collocation points {𝑥𝑖}𝑁𝑖=0 and 𝐹𝑚(𝑦) is the 𝑚th Lagrange interpolation polynomial associated with the collocation points {𝑦𝑖}𝑀𝑖=0.

Using a (𝑀+1)-point Gauss quadrature formula relative to the Jacobi weights {𝜔𝑗}𝑀𝑗=0, (3.12) can be approximated by 𝑣𝑗𝑖𝑥=𝑔𝑖,𝑦𝑗+𝑇1+𝑦𝑗/2𝜇22𝜇1Γ(𝜇)𝑣𝑀𝑁𝑥𝑖,𝑠𝑗(),1𝑀𝑥=𝑔𝑖,𝑦𝑗+𝑇1+𝑦𝑗/2𝜇22𝜇1×Γ(𝜇)𝑀𝑘=0𝑁𝑀𝑛=0𝑚=0𝑣𝑚𝑛𝜆2𝑛𝑥𝑖𝜈𝑛𝑥𝑖+𝐷𝑛𝑥𝑖𝐹𝑚𝑠𝑗𝜃𝑘𝜔𝑘,𝑖=1,,𝑁1;𝑗=0,1,,𝑀,(3.14) where the set {𝜃𝑘}𝑀𝑘=0 coincides with the collocation points {𝑦𝑗}𝑀𝑗=0 on [1,1].

Then the collocation spectral method is to seek 𝑣𝑀𝑁(𝑥,𝑦) of the form (3.13) such that 𝑣𝑗𝑖 satisfies the above collocation equations (3.14) for 1𝑖𝑁1,0𝑗𝑀.

4. Numerical Results with a Collocation Spectral Approximation

In order to demonstrate the effectiveness of the proposed time-space collocation spectral method, some examples are now presented with Dirichlet boundary conditions, Neumann boundary conditions, and mixed boundary conditions.

For completeness sake, the implementation is briefly described here. To simplify the computation, we rewrite the above collocation equations (3.14) into the following: 𝑣𝑗𝑖=𝑔𝑗𝑖+𝑁𝑀𝑛=0𝑚=0𝑣𝑚𝑛𝜆2𝑛𝑥𝑖𝜈𝑛𝑥𝑖+𝐷𝑛𝑥𝑖𝑏𝑗𝑀𝑘=0𝐹𝑚𝑠𝑗𝜃𝑘𝜔𝑘=𝑔𝑗𝑖+𝑁𝑀𝑛=0𝑚=0𝑑𝑖𝑛𝑎𝑗𝑚𝑣𝑚𝑛=𝑔𝑗𝑖+𝑁1𝑀𝑛=1𝑚=0𝑑𝑖𝑛𝑎𝑗𝑚𝑣𝑚𝑛+𝑀𝑚=0𝑑𝑖0𝑣𝑚0+𝑑𝑖𝑁𝑣𝑚𝑁𝑎𝑗𝑚,𝑖=1,,𝑁1;𝑗=0,1,,𝑀,(4.1) where 𝑏𝑗=𝑇1+𝑦𝑗/2𝜇22𝜇1Γ(𝜇),𝑎𝑗𝑚=𝑏𝑗𝑀𝑘=0𝐹𝑚𝑠𝑗𝜃𝑘𝜔𝑘,𝑔𝑗𝑖𝑥=𝑔𝑖,𝑦𝑗,𝑑𝑖𝑛=𝜆2𝑛𝑥𝑖𝜈𝑛𝑥𝑖+𝐷𝑛𝑥𝑖.(4.2)

In our numerical tests, we use the Chebyshev Gauss-Lobatto collocation points 𝑥𝑖1=𝑎+21cos𝑖𝜋𝑁(𝑏𝑎),𝑖=0,1,,𝑁,(4.3) with the associated weights 𝜔0=𝜔𝑁=(𝑏𝑎)𝜋,4𝑁𝜔𝑖=(𝑏𝑎)𝑖𝜋2𝑁,𝑖=1,,𝑁1,(4.4) in the space. The other kinds Gauss-Lobatto collocation points (such as Legendre Gauss-Lobatto collocation points) also can be used. The advantage of Gauss-Lobatto points is that they include the boundary points, which means we can apply boundary conditions there. For the time, the Jacobi Gauss collocation points are used for {𝑦𝑗}𝑀𝑗=0 with the associated weights {𝜔𝑗}𝑀𝑗=0 and other kinds of the Jacobi collocation points also suit to be used.

Let us set 𝐕𝐌𝐍=𝑣01,𝑣02,,𝑣0𝑁1,𝑣11,,𝑣1𝑁1,,𝑣𝑀𝑁1𝑇;𝐆𝐌𝐍=𝑔01,𝑔02,,𝑔0𝑁1,𝑔11,,𝑔1𝑁1,,𝑔𝑀𝑁1𝑇;(4.5) let 𝐀=(𝑎𝑖𝑗) be a matrix of (𝑀+1) by (𝑀+1), and 𝐃=(𝑑𝑖𝑗) is a matrix of (𝑁1) by (𝑁1).

4.1. Implementation of Dirichlet Boundary Conditions

The Dirichlet boundary conditions are directly applied in (4.1) and give numerical solutions on boundary in the following way: 𝑣𝑚0=𝑇1+𝑦𝑚21𝜇𝜑1𝑇1+𝑦𝑚2,𝑣𝜙(𝑎)𝑚𝑁=𝑇1+𝑦𝑚21𝜇𝜑2𝑇1+𝑦𝑚2.𝜙(𝑏)(4.6)

We set 𝐆𝐌𝐍=𝑔01,𝑔02,,𝑔0𝑁1,𝑔11,,𝑔1𝑁1,,𝑔𝑀𝑁1𝑇,(4.7) where 𝑔𝑗𝑖=𝑀𝑚=0𝑑𝑖0𝑣𝑚0+𝑑𝑖𝑁𝑣𝑚𝑁𝑎𝑗𝑚.(4.8)

Thus, the numerical scheme (4.1) leads to a system of equation of the form 𝐕𝐌𝐍=𝐅𝐌𝐍+𝐂𝐕𝐌𝐍,(4.9) where 𝐅𝐌𝐍=𝐆𝐌𝐍+𝐆𝐌𝐍, 𝐂=[𝑎𝑖𝑗𝐃] is a matrix of (𝑁1)×(𝑀+1) by (𝑁1)×(𝑀+1).

4.2. Implementation of Neumann Boundary Conditions

The Neumann boundary conditions (2.3) at 𝑥=𝑎 and 𝑥=𝑏 can be approximated as 𝜕𝑥𝑣𝑀𝑁𝑥0,𝑦𝑚=𝑁𝑛=0𝑣𝑚𝑛𝑛𝑥0=𝑇1+𝑦𝑚21𝜇𝜑1𝑇1+𝑦𝑚2𝜙𝑥0𝑟𝑚(1),𝜕𝑥𝑣𝑀𝑁𝑥𝑁,𝑦𝑚=𝑁𝑛=0𝑣𝑚𝑛𝑛𝑥𝑁=𝑇1+𝑦𝑚21𝜇𝜑2𝑇1+𝑦𝑚2𝜙𝑥𝑁𝑟𝑚(2).(4.10) Equation (4.10) can be written as follows: 0𝑥0𝑣𝑚0+𝑁𝑥0𝑣𝑚𝑁=𝑟𝑚(1)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥𝑛,0𝑥𝑁𝑣𝑚0+𝑁𝑥𝑁𝑣𝑚𝑁=𝑟𝑚(2)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥𝑛.(4.11) Solving (4.11) for 𝑣𝑚0 and 𝑣𝑚𝑁, we get 𝑣𝑚0=𝑁𝑥0𝑟𝑚(2)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥𝑁𝑁𝑥𝑁𝑟𝑚(1)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥0𝑁𝑥00𝑥𝑁0𝑥0𝑁𝑥𝑁,𝑣𝑚𝑁=0𝑥𝑁𝑟𝑚(1)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥00𝑥0𝑟𝑚(2)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥𝑁𝑁𝑥00𝑥𝑁0𝑥0𝑁𝑥𝑁.(4.12)

Then the last right term of (4.1) can be written as follows: 𝑀𝑚=0𝑑𝑖0𝑣𝑚0+𝑑𝑖𝑁𝑣𝑚𝑁𝑎𝑗𝑚=𝑀𝑚=0𝑐𝑖(1)𝑟𝑚(1)𝑐𝑖(2)𝑟𝑚(2)𝑎𝑗𝑚+𝑀𝑚=0𝑐𝑖(2)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥𝑁𝑐𝑖(1)𝑁1𝑛=1𝑣𝑚𝑛𝑛𝑥0𝑎𝑗𝑚=𝑀𝑚=0𝑐𝑖(1)𝑟𝑚(1)𝑐𝑖(2)𝑟𝑚(2)𝑎𝑗𝑚+𝑀𝑚=0𝑁1𝑛=1𝑣𝑚𝑛𝑐𝑖(2)𝑛𝑥𝑁𝑐𝑖(1)𝑛𝑥0𝑎𝑗𝑚,(4.13) where 𝑐𝑖(1)=0𝑥𝑁𝑑𝑖𝑁𝑁𝑥𝑁𝑑𝑖0𝑁𝑥00𝑥𝑁0𝑥0𝑁𝑥𝑁,𝑐𝑖(2)=0𝑥0𝑑𝑖𝑁𝑁𝑥0𝑑𝑖0𝑁𝑥00𝑥𝑁0𝑥0𝑁𝑥𝑁.(4.14)

Let 𝐆𝐌𝐍=̃𝑔01,̃𝑔02,,̃𝑔0𝑁1,̃𝑔11,,̃𝑔1𝑁1,,̃𝑔𝑀𝑁1𝑇,(4.15) where ̃𝑔𝑗𝑖=𝑀𝑚=0𝑐𝑖(1)𝑟𝑚(1)𝑐𝑖(2)𝑟𝑚(2)𝑎𝑗𝑚,(4.16)𝑑𝑖𝑛=𝑑𝑖𝑛+𝑐𝑖(2)𝑛𝑥𝑁𝑐𝑖(1)𝑛𝑥0=𝜆2𝑛𝑥𝑖𝜈𝑛𝑥𝑖+𝐷𝑛𝑥𝑖+𝑐𝑖(2)𝑛𝑥𝑁𝑐𝑖(1)𝑛𝑥0,𝑑𝐃=𝑖𝑗isamatrixof(𝑁1)by(𝑁1).(4.17)

Then the numerical scheme (4.1) leads to a system of equation of the form 𝐕𝐌𝐍=𝐅𝐌𝐍+𝐂𝐕𝐌𝐍,(4.18) where 𝐅𝐌𝐍=𝐆𝐌𝐍+𝐆𝐌𝐍, 𝐂=[𝑎𝑖𝑗𝐃] is a matrix of (𝑁1)×(𝑀+1) by (𝑁1)×(𝑀+1).

Remark 4.1. The implementation for the mixed boundary conditions also can be derived by (4.6) and (4.12).

4.3. Numerical Experiments

In this subsection, the proposed numerical scheme is applied to several test problems to show the efficiency and spectral accuracy.

In each example, we have calculated 𝐿2 errors and 𝐿 errors given by the following formulas: 𝐿2error=𝑁𝑀𝑖=0𝑗=0𝑢𝑗𝑖𝑥𝑢𝑖,𝑡𝑗2𝑤𝑖𝑤𝑗;𝐿error=max𝑖,𝑗||𝑢𝑗𝑖𝑥𝑢𝑖,𝑡𝑗||,(4.19) where 𝑢𝑗𝑖 is the numerical approximation solutions of the exact solutions 𝑢(𝑥𝑖,𝑡𝑗).

Example 4.2 (Dirichlet Boundary Conditions). In this example, we consider the following time fractional diffusion equations: 𝐷𝜇𝑡𝜕𝑢(𝑥,𝑡)=2𝑢(𝑥,𝑡)𝜕𝑥2+𝑓(𝑥,𝑡),𝑥(1,1),𝑡(0,𝑇)𝑢(𝑥,0)=0,𝑥(𝑎,𝑏)𝑢(1,𝑡)=𝑢(1,𝑡)=0,𝑡(0,𝑇),(4.20) where 𝑓(𝑥,𝑡)=Γ(1+𝛽)𝑡Γ(1+𝛽𝜇)𝛽𝜇sin(2𝜋𝑥)+4𝜋2𝑡𝛽sin(2𝜋𝑥).(4.21) The exact solution is given by 𝑢(𝑥,𝑡)=𝑡𝛽sin(2𝜋𝑥),𝛽>1.(4.22)

In Figure 1, we plot the exact and numerical solutions, and, in Figure 2, we represent the associated errors field. Here we set 𝑁=𝑀=22, 𝜇=0.8, and 𝛽=0.2. The results in Figures 1 and 2 denote that the numerical solution using the proposed collocation spectral method is excellent in agreement with the exact solution at the whole domain. The efficiency of this collocation spectral method can be further confirmed by Figures 34, which are the comparison of the exact solution and numerical solution when the space variable 𝑥 or time variable 𝑡 is fixed for various 𝜇 and 𝛽; here 𝑁=𝑀=22.

fig1
Figure 1: Exact and numerical solution with, 𝜇=0.8,𝛽=0.2.
495202.fig.002
Figure 2: Error filed for 𝜇=0.8,𝛽=0.2.
fig3
Figure 3: (a) The comparison of the exact solution and numerical solution when 𝑡=3. (b) The comparison of the exact solution and numerical solution when 𝑥=0.7557.
fig4
Figure 4: (a) The comparison of the exact solution and numerical solution when 𝑡=3. (b) The compare of the exact solution and numerical solution when 𝑥=0.7557.

The main purpose of the numerical test is to check the convergence behavior of numerical solutions with respect to the polynomial degrees 𝑀 and 𝑁 for several 𝜇, especially the convergence in time because of the fractional derivative in time. In order to investigate the spatial accuracy when 𝑁 increases, we take 𝑀 larger enough so that the time discretization errors are negligible compared with the spatial discretization errors. Similary, for the temporal accuracy, we must keep 𝑁 large enough to preclude spatial errors. 𝐿2 errors and 𝐿 errors in semilog scale varying with the polynomial degree 𝑁 (see (a)) or 𝑀 (see (b)) are shown in Figures 5 and 6 for 𝜇=0.05,0.75 with 𝛽=2.2, and Figure 7 for 𝜇=1 with 𝛽=1.2. It is clear that the spectral convergence is achieved both of spatial and temporal errors. This indicates that the convergence in space and time of the time-space collocation spectral method is exponential, even though for the classical diffusion equation (𝜇=1).

fig5
Figure 5: y-axis is log scale. (a) Errors versus 𝑁 with 𝑀=16,𝜇=0.05. (b) Errors versus 𝑀 with 𝑁=16,𝜇=0.05.
fig6
Figure 6: y-axis is log scale. (a) Errors versus 𝑁 with 𝑀=16,𝜇=0.75. (b) Errors versus 𝑀 with 𝑁=16,𝜇=0.75.
fig7
Figure 7: y-axis is log scale. (a) Errors versus 𝑁 with 𝑀=16,𝜇=1. (b) Errors versus 𝑀 with 𝑁=16,𝜇=1.

Example 4.3 (Neumann Boundary Conditions). Consider the time fractional advection-diffusion equation with Neumann boundary conditions 𝐷𝜇𝑡𝑢(𝑥,𝑡)=𝜈𝜕𝑢(𝑥,𝑡)𝜕𝜕𝑥+𝐷2𝑢(𝑥,𝑡)𝜕𝑥2𝜕+𝑓(𝑥,𝑡),𝑥(3,5),𝑡>0,𝑢(𝑥,0)=0,𝑥(3,5),𝑥𝑢(3,𝑡)=𝜋𝑡22,𝜕𝑥𝑢(5,𝑡)=𝜋𝑡22,𝑡>0.(4.23) We choose the suitable source term 𝑓(𝑥,𝑡) to obtain the exact solution 𝑢(𝑥,𝑡)=𝑡2cos𝜋𝑥2.(4.24)

In Table 1, 𝐿2 errors and 𝐿 errors varying with (𝑁,𝑀) and 𝜇 are given, and it indicates the efficiency of the present technique.

tab1
Table 1: 𝐿2 errors and 𝐿 errors varying with (𝑁, 𝑀) and 𝜇.

In Figure 8, we plot 𝐿2 errors and 𝐿 errors in semilog scale. Figures 8(a) and 8(b) are concerned with the spatial errors with 𝑀=16 and the temporal errors with 𝑁=16, respectively. As expected, the spatial and temporal spectral convergence is achieved. Our greatest interest is to check the convergence in time because of the fractional derivative in time. So we plot the 𝐿2 errors and 𝐿 errors as functions of 𝑀 for several more values 𝜇 in Figure 9, where 𝑁=16. This graph shows that the spectral convergence in time can be reached, and that the time-space collocation spectral method also works well for the classical advection-diffusion equation.

fig8
Figure 8: y-axis is log scale. (a) Errors versus 𝑁 with 𝑀=16,𝜇=0.9. (b) Errors versus 𝑀 with 𝑁=16,𝜇=0.9.
fig9
Figure 9: y-axis is log scale. Errors versus 𝑀 with 𝑁=16 for several 𝜇.

Example 4.4 (Mixed Boundary Conditions). Consider the time fractional reaction-subdiffusion equation 𝐷𝜇𝑡3𝑢(𝑥,𝑡)=4𝜕𝑢(𝑥,𝑡)+2𝑢(𝑥,𝑡)𝜕𝑥2𝑥,𝑥(0,2𝜋),𝑡>0,𝑢(𝑥,0)=sin2,𝑥(0,2𝜋);𝑢(0,𝑡)=0,𝜕𝑥1𝑢(2𝜋,𝑡)=2𝐸𝜇(𝑡𝜇).(4.25) The exact solution is given by 𝑢(𝑥,𝑡)=𝐸𝜇(𝑡𝜇𝑥)sin2.(4.26)

In Figure 10, we plot the visual fields of the exact and numerical solutions. Comparing Figures 10(a) and 10(c), and Figures 10(b) and 10(d), we see that the numerical solution is in good agreement with the exact solution; the shapes on the space variable are similar to Sine function. In addition, we can see that the solutions decay over time from the initial state. Comparing Figures 10(a) and 10(b) with Figures 10(c) and 10(d), (a) and (b) decay faster in the beginning and the trend is the opposite of (c) and (d). This is consistent with the behavior of the function 𝐸𝛼(𝑧). For 0<𝑡<0.77, 𝐸1/2(𝑡1/2) decays faster than 𝐸1(𝑡), whereas for 𝑡>0.77, the trend is just the opposite.

fig10
Figure 10: Exact and numerical solutions.

We plot in Figure 11 the 𝐿2 errors and 𝐿 errors versus 𝑁 with 𝑀=16 in 11(a) and versus 𝑀 with 𝑁=16 in 11(b). As can be seen in Figure 11, the proposed method provides spatial and temporal spectral convergence for 𝐿2 errors and 𝐿 errors.

fig11
Figure 11: y-axis is log scale. (a) Error varying with the polynomial degree 𝑁. (b) Error varying with the polynomial degree 𝑀.

5. Conclusion

This paper proposes a time-space collocation spectral method for a class of time fractional differential equations with Caputo derivatives. The proposed method also works well for the corresponding classical integer-order partial differential equations (𝜇=1), and it differs from (and is simpler than) the existing time-space spectral methods which are based on the Petrov-Galerkin or Dual-Petrov-Galerkin formulation. The main advantage of the present scheme is that it gives very accurate convergency by choosing less number of grid points and the problem can be solved up to big time, and the storage requirement due to the time memory effect can be considerably reduced. At the same time, the technique is also simpler and easier to apply to multidimensional problems than the existing Galerkin spectral method such as the methods in [19, 20].

References

  1. I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1999. View at Zentralblatt MATH
  2. W. Wyss, “The fractional diffusion equation,” Journal of Mathematical Physics, vol. 27, no. 11, pp. 2782–2785, 1986. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  3. W. R. Schneider and W. Wyss, “Fractional diffusion and wave equations,” Journal of Mathematical Physics, vol. 30, no. 1, pp. 134–144, 1989. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  4. F. Huang and F. Liu, “The time fractional diffusion equation and the advection-dispersion equation,” The Australian & New Zealand Industrial and Applied Mathematics Journal, vol. 46, no. 3, pp. 317–330, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. F. Huang and B. Guo, “General solutions to a class of time fractional partial differential equations,” Applied Mathematics and Mechanics, vol. 31, no. 7, pp. 815–826, 2010.
  6. Y. Luchko and R. Goren°o, “The initial value problem for some fractional differen-tial equations with the Caputo derivatives,” Preprint Series A08-98, Fachbreich Mathematik und Informatik, Freic Universitat Berlin, 1998.
  7. N. T. Shawagfeh, “Analytical approximate solutions for nonlinear fractional differential equations,” Applied Mathematics and Computation, vol. 131, no. 2-3, pp. 517–529, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  8. S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives, Gordon and Breach Science Publishers, Yverdon, Switzerland, 1993.
  9. J. Chen, F. Liu, and V. Anh, “Analytical solution for the time-fractional telegraph equation by the method of separating variables,” Journal of Mathematical Analysis and Applications, vol. 338, no. 2, pp. 1364–1377, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for fractional advection-dispersion flow equations,” Journal of Computational and Applied Mathematics, vol. 118, pp. 283–299, 2004.
  11. T. A. M. Langlands and B. I. Henry, “The accuracy and stability of an implicit solution method for the fractional diffusion equation,” Journal of Computational Physics, vol. 205, no. 2, pp. 719–736, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. Y. Lin and C. Xu, “Finite difference/spectral approximations for the time-fractional diffusion equation,” Journal of Computational Physics, vol. 225, no. 2, pp. 1533–1552, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. C.-M. Chen, F. Liu, and K. Burrage, “Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation,” Applied Mathematics and Computation, vol. 198, no. 2, pp. 754–769, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  14. I. Podlubny, A. Chechkin, T. Skovranek, Y. Chen, and B. M. Vinagre Jara, “Matrix approach to discrete fractional calculus. II. Partial fractional differential equations,” Journal of Computational Physics, vol. 228, no. 8, pp. 3137–3153, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  15. G. J. Fix and J. P. Roop, “Least squares finite-element solution of a fractional order two-point boundary value problem,” Computers & Mathematics with Applications, vol. 48, no. 7-8, pp. 1017–1033, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. J. P. Roop, “Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in R2,” Journal of Computational and Applied Mathematics, vol. 193, no. 1, pp. 243–268, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  17. W. H. Deng, “Finite element method for the space and time fractional Fokker-Planck equation,” SIAM Journal on Numerical Analysis, vol. 47, no. 1, pp. 204–226, 2008/09. View at Publisher · View at Google Scholar
  18. E. Hanert, “A comparison of three Eulerian numerical methods for fractional-order transport models,” Environmental Fluid Mechanics, vol. 10, no. 1, pp. 7–20, 2010. View at Publisher · View at Google Scholar · View at Scopus
  19. E. Hanert, “On the numerical solution of space-time fractional diffusion models,” Computers and Fluids, vol. 46, no. 1, pp. 33–39, 2011. View at Publisher · View at Google Scholar · View at Scopus
  20. X. Li and C. Xu, “A space-time spectral method for the time fractional diffusion equation,” SIAM Journal on Numerical Analysis, vol. 47, no. 3, pp. 2108–2131, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. Y. Chen and T. Tang, “Convergence analysis of the Jacobi spectral-collocation methods for Volterra integral equations with a weakly singular kernel,” Mathematics of Computation, vol. 79, no. 269, pp. 147–167, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH