Abstract

We consider finite element Galerkin solutions for the space fractional diffusion equation with a nonlinear source term. Existence, stability, and order of convergence of approximate solutions for the backward Euler fully discrete scheme have been discussed as well as for the semidiscrete scheme. The analytical convergent orders are obtained as 𝑂(𝑘+̃𝛾), where ̃𝛾 is a constant depending on the order of fractional derivative. Numerical computations are presented, which confirm the theoretical results when the equation has a linear source term. When the equation has a nonlinear source term, numerical results show that the diffusivity depends on the order of fractional derivative as we expect.

1. Introduction

Fractional calculus is an old mathematical topic but it has not been attracted enough for almost three hundred years. However, it has been recently proven that fractional calculus is a significant tool in the modeling of many phenomena in various fields such as engineering, physics, porous media, economics, and biological sciences. One can see related references in [17].

In the classical diffusion model, it is assumed that particles are distributed in a normal bell-shaped pattern based on the Brownian motion. In general, the nature of diffusion is characterized by the mean squared displacement (𝑟)2=2𝑑𝜅𝜇𝑡𝜇,(1.1) where 𝑑 is the spatial dimension and 𝜅𝜇 is the diffusion constant. The classical normal diffusion case arises when the exponent 𝜇=1. When 𝜇1, anomalous diffusions arise. The anomalous diffusion is classified as the process is subdiffusive (diffusive slowly) when 𝜇<1 or superdiffusive (diffusive fast) when 𝜇>1.

As mentioned before, in many real problems, it is more adequate to use anomalous diffusion described by fractional derivatives than the classical normal diffusion [4, 5, 812]. One typical model for anomalous diffusion is the fractional superdiffusion equation arising in chaotic and turbulent processes, where the usual second derivative in space is replaced by a fractional derivative of order 1<𝜇<2.

In this paper we discuss Galerkin approximate solutions for the space fractional diffusion equation with a nonlinear source term. The equation is described as 𝜕𝑢(𝑥,𝑡)𝜕𝑡=𝜅𝜇𝜇𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡,𝑢)(1.2) with an initial condition 𝑢(𝑥,0)=𝑢0,𝑥Ω𝐑(1.3) and boundary conditions 𝑢(𝑥,𝑡)=0,𝑥𝜕Ω,0𝑡𝑇,(1.4) where 𝜅𝜇 denotes the anomalous diffusion coefficient and 𝜕Ω is the boundary of the domain Ω. And the differential operator 𝜇 is 𝜇=12𝑎𝐷𝜇𝑥+12𝑥𝐷𝜇𝑏,(1.5) where 𝑎𝐷𝜇𝑥 and 𝑥𝐷𝜇𝑏 are called the left and the right Riemann-Liouville space fractional derivatives of order 𝜇, respectively, defined by 𝐃𝜇𝑢=𝑎𝐷𝜇𝑥𝑢(𝑥)=𝐷𝑛𝑎𝐷𝑥𝜇𝑛1𝑢(𝑥)=𝑑Γ(𝑛𝜇)𝑛𝑑𝑥𝑛𝑥𝑎(𝑥𝜉)𝑛𝜇1𝐃𝑢(𝜉)𝑑𝜉,𝜇𝑢=𝑥𝐷𝜇𝑏𝑢(𝑥)=(𝐷)𝑛𝑥𝐷𝑏𝜇𝑛𝑢(𝑥)=(1)𝑛Γ𝑑(𝑛𝜇)𝑛𝑑𝑥𝑛𝑏𝑥(𝜉𝑥)𝑛𝜇1𝑢(𝜉)𝑑𝜉.(1.6) Here 𝑛 is the smallest integer such that 𝑛1𝜇<𝑛.

Throughout this paper, we will assume that the nonlinear source term 𝑓(𝑥,𝑡,𝑢) is locally Lipschitz continuous with constants 𝐶𝑙 and 𝐶𝑓 such that 𝑓(𝑢)𝑓(𝑣)𝐿2(Ω)𝐶𝑙𝑢𝑣𝐿2(Ω),(1.7)𝑓(𝑢)𝐿2(Ω)𝐶𝑓𝑢𝐿2(Ω)(1.8) for 𝑢,𝑣{𝑤𝐻0𝜇/2(Ω)𝑤𝐿2(Ω)𝑙}.

Baeumer et al. [8, 13] have proved existence and uniqueness of a strong solution for (1.2) using the semigroup theory when 𝑓(𝑥,𝑡,𝑢) is globally Lipschitz continuous. Furthermore, when 𝑓(𝑥,𝑡,𝑢) is locally Lipschitz continuous, existence of a unique strong solution has also been shown by introducing the cut-off function.

Finite difference methods have been studied in [1416] for linear space fractional diffusion problems. They used the right-shifted Grüwald-Letnikov approximate for the fractional derivative since the standard Grüwald-Letnikov approximate gives the unconditional instability even for the implicit method. Using the right-shifted Grüwald-Letnikov approximation, the method of lines has been applied in [12] for numerical approximate solutions.

For the space fractional diffusion problems with a nonlinear source term, Lynch et al. [17] used the so-called L2 and L2C methods in [6] and compared computational accuracy of them. Baeumer et al. [8] give existence of the solution and computational results using finite difference methods. Choi et al. [18] have shown existence and stability of numerical solutions of an implicit finite difference equation obtained by using the right-shifted Grüwald-Letnikov approximation. For the time fractional diffusion equations, explicit and implicit finite difference methods have been used in [11, 1923].

Compared to finite difference methods on the fractional diffusion equation, finite element methods have been rarely discussed. Ervin and Roop [24] have considered finite element analysis for stationary linear advection dispersion equations, and Roop [25] has studied finite element analysis for nonstationary linear advection dispersion equations. The finite element numerical approximations have been discussed for the time and space fractional Fokker-Planck equation in Deng [9] and for the space general fractional diffusion equations with a nonlocal quadratic nonlinearity but a linear source term in Ervin et al. [26].

As far as we know, finite element methods have not been considered for the space fractional diffusion equation with nonlinear source terms. In this paper, we will discuss finite element solutions for the problem (1.2)–(1.4) under the assumption of existence of a sufficiently regular solution 𝑢 of the equation. Finite element numerical analysis of the semidiscrete and fully discrete methods for (1.2)–(1.4) will be considered using the backward Euler method in time and Galerkin finite element method in space as well as the semidiscrete method. We will discuss existence, uniqueness, and stability of the numerical solutions for the problem (1.2)–(1.4). Also, 𝐿2-error estimate will be considered for the problem (1.2)–(1.4).

The outline of the paper is as follows. We introduce some properties of the space fractional derivatives in Section 2, which will be used in later discussion. In Section 3, the semidiscrete variational formulation for (1.2) based on Galerkin method is given. Existence, stability and 𝐿2-error estimate of the semidiscrete solution are analyzed. In Section 4, existence and unconditional stability of approximate solutions for the fully discrete backward Euler method are shown following the idea of the semidiscrete method. Further, 𝐿2-error estimates are obtained, whose convergence is of )𝑂(𝑘+̃𝛾, where ̃𝛾=𝜇 if 𝜇3/2 and ̃𝛾=𝜇𝜖,0<𝜖<1/2, if 𝜇=3/2. Finally, numerical examples are given in order to see the theoretical convergence order discussed in Section 5. We will see that numerical solutions of fractional diffusion equations diffuse more slowly than that of the classical diffusion problem and diffusivity depends on the order of fractional derivatives.

2. The Variational Form

In this section we will consider the variational form of problem (1.2)–(1.4) and show existence and stability of the weak solution. We first recall some basic properties of Riemann-Liouville fractional calculus [9, 24].

For any given positive number 𝜇>0, define the seminorm |𝑢|𝐽𝜇𝐿(𝐑)=𝐃𝜇𝑢𝐿2(𝐑)(2.1) and the norm 𝑢𝐽𝜇𝐿(𝐑)=𝑢2𝐿2(𝐑)+|𝑢|2𝐽𝜇𝐿(𝐑)1/2,(2.2) where the left fractional derivative space 𝐽𝜇𝐿(𝐑) denotes the closure of 𝐶0(𝐑) with respect to the norm 𝐽𝜇𝐿(𝐑).

Similarly, we may define the right fractional derivative space 𝐽𝜇𝑅(𝐑) as the closure of 𝐶0(𝐑) with respect to the norm 𝐽𝜇𝑅(𝐑), where 𝑢𝐽𝜇𝑅(𝐑)=𝑢2𝐿2(𝐑)+|𝑢|2𝐽𝜇𝑅(𝐑)1/2(2.3) and the seminorm |𝑢|𝐽𝜇𝑅(𝐑)=𝐃𝜇𝑢𝐿2(𝐑).(2.4)

Furthermore, with the help of Fourier transform we define a seminorm |𝑢|𝐻𝜇(𝐑)=|𝜔|𝜇̂𝑢𝐿2(𝐑)(2.5) and the norm 𝑢𝐻𝜇(𝐑)=𝑢2𝐿2(𝐑)+|𝑢|2𝐻𝜇(𝐑)1/2.(2.6) Here 𝐻𝜇(𝐑) denotes the closure of 𝐶0(𝐑) with respect to 𝐻𝜇(𝐑). It is known in [24] that the spaces 𝐽𝜇𝐿(𝐑),𝐽𝜇𝑅(𝐑), and 𝐻𝜇(𝐑) are all equal with equivalent seminorms and norms. Analogously, when the domain Ω is a bounded interval, the spaces 𝐽𝜇𝐿,0(Ω),𝐽𝜇𝑅,0(Ω), and 𝐻𝜇0(Ω) are equal with equivalent seminorms and norms [24, 27].

The following lemma on the Riemann-Liouville fractional integral operators will be used in our analysis, which can be proved by using the property of Fourier transform [24].

Lemma 2.1. For a given 𝜇>0 and a real valued function 𝑢𝐃𝜇𝑢,𝐃𝜇𝑢=cos(𝜋𝜇)𝐃𝜇𝑢2𝐿2(𝐑).(2.7)

Remark 2.2. It follows from (2.7) that we may use the following norm: 𝑢2𝐻0𝜇/2(𝐑)=𝑢2𝐿2(𝐑)+𝜅𝜇|||𝜇cos𝜋2||||𝑢|2𝐻0𝜇/2(𝐑)(2.8) instead of the norm 𝑢𝐻𝜇(𝐑).

For the seminorm on 𝐻𝜇0(Ω) with Ω=(𝑎,𝑏), the following fractional Poincaré-Friedrich’s inequality holds. For the proof, we refer to [9, 24].

Lemma 2.3. For 𝑢𝐻𝜇0(Ω), there is a positive constant 𝐶 such that 𝑢𝐿2(Ω)𝐶|𝑢|𝐻𝜇0(Ω)(2.9) and for 0<𝑠<𝜇,𝑠𝑛1/2,𝑛1𝜇<𝑛,𝑛, |𝑢|𝐻𝑠0(Ω)𝐶|𝑢|𝐻𝜇0(Ω).(2.10)

Hereafter, a positive number 𝐶 will denote a generic constant. Also the semigroup property and the adjoint property hold for the Riemann-Liouville fractional integral operators [9, 24]: for all 𝜇,𝜈>0, if 𝑢𝐿𝑝(Ω), 𝑝1, then 𝑎𝐷𝑥𝑎𝜇𝐷𝑥𝜈𝑢(𝑥)=𝑎𝐷𝑥𝜇𝜈𝑢(𝑥),𝑥Ω,𝑥𝐷𝑏𝑥𝜇𝐷𝑏𝜈𝑢(𝑥)=𝑥𝐷𝑏𝜇𝜈𝑢(𝑥),𝑥Ω,(2.11) and specially 𝑎𝐷𝑥𝜇𝑢,𝑣𝐿2(Ω)=𝑢,𝑥𝐷𝑏𝜇𝑣𝐿2(Ω),𝑢,𝑣𝐿2(Ω).(2.12)

In the rest of this section, we will consider a weak problem for (1.2)–(1.4) with 1<𝜇<2: find a function 𝑢𝐻0𝜇/2(Ω) such that 𝑢𝑡=𝜅,𝑣𝜇𝜇𝑢,𝑣+(𝑓(𝑢),𝑣),𝑣𝐻0𝜇/2(Ω).(2.13)

Since there is a weak solution of (2.13) when 𝑓 is locally Lipschitz continuous as in [8, 13], we here only discuss the stability of the weak solution, to show that we need the following lemma.

Lemma 2.4. For all 𝑣𝐻0𝜇/2(Ω), the following inequality holds: 𝜅𝜇𝜇𝑣,𝑣𝜅𝜇|||𝜇cos𝜋2||||𝑣|2𝐻0𝜇/2(Ω).(2.14)

Proof. Following the ideas in [9, 26], we obtain the following inequality by using the properties (2.11)-(2.12) and Lemmas 2.12.3: 𝜅𝜇𝜇𝜅𝑣,𝑣=𝜇2𝑎𝐷𝜇𝑥+𝑣,𝑣𝑥𝐷𝜇𝑏𝜅𝑣,𝑣=𝜇2𝑏𝑎𝐷2𝑎𝐷𝑥(2𝜇)𝑣𝑣𝑑𝑥+𝑏𝑎(𝐷)2𝑥𝐷𝑏(2𝜇)𝑣=𝜅𝑣𝑑𝑥𝜇2𝑏𝑎𝐷𝑎𝐷𝑥(2𝜇)𝑣𝐷𝑣𝑑𝑥+𝑏𝑎𝐷𝑥𝐷𝑏(2𝜇)𝑣=𝜅𝐷𝑣𝑑𝑥𝜇2𝑏𝑎𝑎𝐷𝑥(2𝜇)𝐷𝑣𝐷𝑣𝑑𝑥+𝑏𝑎𝑥𝐷𝑏(2𝜇)=𝜅𝐷𝑣𝐷𝑣𝑑𝑥𝜇2𝑏𝑎𝑎𝐷𝑥𝑎(2𝜇)/2𝐷𝑥(2𝜇)/2+𝐷𝑣𝐷𝑣𝑑𝑥𝑏𝑎𝑥𝐷𝑏𝑥((2𝜇)/2)𝐷𝑏(2𝜇)/2=𝜅𝐷𝑣𝐷𝑣𝑑𝑥𝜇2𝑏𝑎𝑎𝐷𝑥(2𝜇)/2𝐷𝑣𝑥𝐷𝑏(2𝜇)/2+𝐷𝑣𝑑𝑥𝑏𝑎𝑥𝐷𝑏(2𝜇)/2𝐷𝑣𝑎𝐷𝑥(2𝜇)/2𝐷𝑣𝑑𝑥=𝜅𝜇𝐃𝜇/2𝑣,𝐃(𝜇/2)𝑣=𝜅𝜇𝜇cos𝜋2𝐃𝜇/2𝑣2𝐿2(Ω)𝜅𝜇|||𝜇cos𝜋2||||𝑣|2𝐻0𝜇/2(Ω).(2.15) This completes the proof.

We consider the stability of a weak solution 𝑢 for (2.13).

Theorem 2.5. Let 𝑢 be a solution of (2.13). Then there is a constant 𝐶 such that 𝑢(𝑡)𝐿2(Ω)𝐶𝑢(0)𝐿2(Ω).(2.16)

Proof. Taking 𝑣=𝑢(𝑡) in (2.13), we obtain 𝑢𝑡𝜅,𝑢𝜇𝜇=𝑢,𝑢(𝑓(𝑢),𝑢).(2.17) Since the second term on the left hand side is nonnegative from Lemma 2.4, we have 12𝑑𝑑𝑡𝑢2𝐿2(Ω)12𝑑𝑑𝑡𝑢2𝐿2(Ω)+𝜅𝜇|||𝜇cos𝜋2||||𝑢|2𝐻0𝜇/2(Ω)𝑓(𝑢)𝐿2(Ω)𝑢𝐿2(Ω)𝐶𝑓𝑢2𝐿2(Ω).(2.18) Integrating both sides with respect to 𝑡, we obtain 𝑢(𝑡)2𝐿2(Ω)𝑢(0)2𝐿2(Ω)+𝐶𝑡0𝑢(𝑠)2𝐿2(Ω)𝑑𝑠.(2.19) An application of Gronwall's inequality gives that there is a constant 𝐶 such that 𝑢(𝑡)2𝐿2(Ω)𝐶𝑢(0)2𝐿2(Ω).(2.20) This completes the proof.

3. The Semidiscrete Variational Form

In this section, we will analyze the stability and error estimates of Galerkin finite element solutions for the semidiscrete variational formulation for (1.2).

Let 𝑆 be a partition of Ω with a grid parameter such that Ω={𝐾𝐾𝑆} and =max𝐾𝑆𝐾, where 𝐾 is the width of the subinterval 𝐾. Associated with the partition 𝑆, we may define a finite-dimensional subspace 𝑉𝐻0𝜇/2(Ω) with a basis {𝜑𝑖}𝑁𝑖=1 of piecewise polynomials. Then the semidiscrete variational problem is to find 𝑢𝑉 such that 𝑢,𝑡=𝜅,𝑣𝜇𝜇𝑢+𝑓𝑢,𝑣,𝑣,𝑣𝑉,𝑢(3.1)(𝑥,0)=𝑢0,𝑢(3.2)(𝑎,𝑡)=𝑢(𝑏,𝑡)=0.(3.3) Since 𝑢 can be represented as 𝑢(𝑥,𝑡)=𝑁𝑖=1𝛼𝑖(𝑡)𝜑𝑖(𝑥),(3.4) we may rewrite (3.1) in a matrix form: 𝐀̇𝐮(𝑡)+𝐁𝐮=𝐅(𝐮),(3.5) where 𝑁×𝑁 matrices 𝐀 and 𝐁 and vectors 𝐮 and 𝐅 are 𝑎𝐀=𝑖𝑗,𝑎𝑖𝑗=𝜑𝑖,𝜑𝑗,𝑏𝐁=𝑖𝑗,𝑏𝑖𝑗𝜅=𝜇2𝐃𝜇/2𝜑𝑖,𝐃(𝜇/2)𝜑𝑗+𝐃𝜇/2𝜑𝑗,𝐃(𝜇/2)𝜑𝑖,𝐹𝐅(𝐮)=𝑗,𝐹𝑗=𝑓𝑁𝑙=1𝛼𝑙𝜑𝑙,𝜑𝑗,𝛼𝐮=1(𝑡),𝛼2(𝑡),,𝛼𝑁(𝑡)𝑇.(3.6)

It follows from 𝑁𝑖,𝑗=1𝛼𝑖𝛼𝑗(𝜑𝑖,𝜑𝑗)=(𝑁𝑖=1𝛼𝑖𝜑𝑖,𝑁𝑗=1𝛼𝑗𝜑𝑗)0 and Lemma 2.4 that matrices 𝐀 and 𝐁 are nonnegative definite and nonsingular. Thus this system (3.5) of ordinary differential equations has a unique solution since 𝑓 is locally Lipschitz continuous.

The stability for the semidiscrete variational problem (3.1) can be obtained by following the proof of Theorem 2.5, which is 𝑢𝐿2(Ω)𝑢𝐶0𝐿2(Ω).(3.7)

Now we will consider estimates of error between the weak solution of (2.13) and the one of semidiscrete form (3.1). The finite dimensional subspace 𝑉𝐻0𝜇/2(Ω) is chosen so that the interpolation 𝐼𝑢 of 𝑢 satisfies an approximation property [9, 28]: for 𝑢𝐻𝛾(Ω), 0<𝛾𝑛, and 0𝑠𝛾, there exists a constant 𝐶 depending only on Ω such that 𝑢𝐼𝑢𝐻𝑠(Ω)𝐶𝛾𝑠𝑢𝐻𝛾(Ω).(3.8) Since the norm 𝐻𝑠(Ω) is equivalent to the seminorm ||𝐻𝑠(Ω), we may replace (3.8) by the relation 𝑢𝐼𝑢𝐻𝑠(Ω)𝐶𝛾𝑠|𝑢|𝐻𝛾(Ω).(3.9)

Further we need an adjoint problem to find 𝑤𝐻𝜇(Ω)𝐻0𝜇/2(Ω) satisfying 𝜅𝜇𝜇𝑤=𝑔,inΩ,𝑤=0,on𝜕Ω.(3.10) Bai and Lü [29] have proved existence of a solution to the problem (3.10). We assume as in Ervin and Roop [24] that the solution 𝑤 satisfies the regularity 𝑤𝐻𝜇(Ω)𝐶𝑔𝐿2(Ω)3,𝜇2,(3.11)𝑤𝐻𝜇𝜖(Ω)𝐶𝑔𝐿2(Ω)3,𝜇=21,0<𝜖<2.(3.12)

Let ̃𝑢=𝑃𝑢 be the elliptic projection 𝑃𝐻0𝜇/2(Ω)𝑉 of the exact solution 𝑢, which is defined by 𝜅𝜇𝜇𝑢̃𝑢,𝑣=0,𝑣𝑉.(3.13) Let 𝜃=𝑢̃𝑢 and 𝜌=̃𝑢𝑢. Then the error is expressed as 𝑒=𝑢𝑢𝑢=̃𝑢+̃𝑢𝑢=𝜃+𝜌.(3.14)

First, we consider the following estimates on 𝜌.

Lemma 3.1. Let ̃𝑢 be a solution of (3.13) and let 𝑢𝐻𝜇(Ω)𝐻0𝜇/2(Ω) be the solution of (2.13). Let 𝜌(𝑡)=̃𝑢(𝑡)𝑢(𝑡). Then there is a constant 𝐶 such that 𝜌(𝑡)𝐿2(Ω)𝐶̃𝛾𝑢(𝑡)𝐻𝛾(Ω),𝜌𝑡(𝑡)𝐿2(Ω)(𝐶̃𝛾𝑢𝑡)𝐻𝛾(Ω),(3.15) where ̃𝛾=𝜇 if 𝜇3/2 and ̃𝛾=𝜇𝜖, 0<𝜖<1/2 if 𝜇=3/2.

Proof. It follows from the fractional Poincaré-Friedrich’s inequality and the adjoint property (2.12) that for 𝜓,𝜒𝑉𝐻0𝜇/2(Ω)(𝐃𝜇𝜓,𝜒)=𝑏𝑎𝐃𝜇/2𝜓𝐃(𝜇/2)||𝜓||𝜒𝑑𝑥𝐽𝜇/2𝐿,0(Ω)||𝜒||𝐽𝜇/2𝑅,0(Ω)𝐶𝜓𝐻0𝜇/2(Ω)𝜒𝐻0𝜇/2(Ω).(3.16) Similarly we obtain 𝐃𝜇=𝜓,𝜒𝑏𝑎𝐃(𝜇/2)𝜓𝐃𝜇/2𝜒𝑑𝑥𝐶𝜓𝐻0𝜇/2(Ω)𝜒𝐻0𝜇/2(Ω).(3.17)
It follows from Lemma 2.4 that for 𝑣𝑉𝜅𝜇|||𝜇cos𝜋2|||||𝑢̃𝑢||2𝐻0𝜇/2(Ω)𝜅𝜇𝜇𝑢̃𝑢,𝑢̃𝑢𝜅𝜇𝜇𝑢̃𝑢,𝑢𝑣𝜅𝜇𝜇𝑢̃𝑢,𝑣̃𝑢𝐶𝑢̃𝑢𝐻0𝜇/2(Ω)𝑢𝑣𝐻0𝜇/2(Ω).(3.18) Using the equivalence of seminorms and norms, we obtain 𝑢̃𝑢𝐻0𝜇/2(Ω)𝐶inf𝑣𝑉𝑢𝑣𝐻0𝜇/2(Ω)𝐶𝑢𝐼𝑢𝐻0𝜇/2(Ω).(3.19)
In case of 𝜇3/2 and 𝑣𝑉, by taking 𝑔=𝜌 in (3.10) and using (3.13), (3.16)-(3.17) and the adjoint property (2.12), we have (𝜌,𝜌)=𝜅𝜇(𝜇𝑤,𝜌)=𝜅𝜇(𝜇(𝑤𝑣),𝜌)𝜅𝜇(𝜇𝜌,𝑣)=𝜅𝜇(𝜇(𝑤𝑣),𝜌)𝐶𝑤𝑣𝐻0𝜇/2(Ω)𝜌𝐻0𝜇/2(Ω).(3.20) Taking 𝑣=𝐼𝑤 in the previously mentioned inequalities, we have 𝜌2𝐿2(Ω)𝐶𝑤𝐼𝑤𝐻0𝜇/2(Ω)𝜌𝐻0𝜇/2(Ω)𝐶𝜇/2𝑤𝐻𝜇(Ω)𝑢𝐼𝑢𝐻0𝜇/2(Ω)𝐶𝜇/2𝜌𝐿2(Ω)𝜇/2𝑢𝐻𝜇(Ω).(3.21) Thus we obtain 𝜌𝐿2(Ω)𝐶𝜇𝑢𝐻𝜇(Ω).(3.22) We now differentiate (3.13). Then we obtain 𝜅𝜇(𝜇𝜌𝑡,𝑣)=0 for all 𝑣𝑉. Using the previous duality arguments again, we have 𝜌𝑡𝐿2(Ω)𝐶𝜇𝑢𝐻𝜇(Ω).(3.23)
In case of 𝜇=3/2, we can similarly prove (3.15) by applying the assumption (3.12). This completes the proof.

We now consider the estimates on 𝜃.

Lemma 3.2. Let 𝑢 and ̃𝑢 be the solutions of (3.1)–(3.3) and (3.13), respectively. Let 𝜃(𝑡)=𝑢(𝑡)̃𝑢(𝑡). Then there is a constant 𝐶 such that 𝜃(𝑡)𝐿2(Ω),𝐶̃𝛾(3.24) where ̃𝛾=𝜇 if 𝜇3/2 and ̃𝛾=𝜇𝜖, 0<𝜖<1/2 if 𝜇=3/2.

Proof. It follows from (3.1) and (3.13) that for 𝑣𝑉, 𝜃𝑡,𝑣𝜅𝜇(𝜇𝑓𝑢𝜃,𝑣)=𝜌𝑓(𝑢),𝑣𝑡.,𝑣(3.25) Replacing 𝑣=𝜃 in (3.25), we obtain 12𝑑𝑑𝑡𝜃2𝐿2(Ω)𝐶𝑙𝑢𝑢𝐿2(Ω)𝜃𝐿2(Ω)+𝜌𝑡𝐿2(Ω)𝜃𝐿2(Ω).(3.26) Using Young's inequality 𝑑𝑑𝑡𝜃2𝐿2(Ω)𝑢𝐶̃𝑢𝐿2(Ω)+̃𝑢𝑢𝐿2(Ω)𝜃𝐿2(Ω)+𝜌𝑡𝐿2(Ω)𝜃𝐿2(Ω)𝐶𝜃𝐿2(Ω)+𝜌𝐿2(Ω)+𝜌𝑡𝐿2(Ω)𝜃𝐿2(Ω)𝐶1𝜃2𝐿2(Ω)+𝐶2𝜌2𝐿2(Ω)+𝐶3𝜌𝑡2𝐿2(Ω).(3.27) Integration on time 𝑡 gives 𝜃(𝑡)2𝐿2(Ω)𝜃(0)2𝐿2(Ω)+𝐶𝑡0𝜃2𝐿2(Ω)𝑑𝑠+𝐶𝑡0𝜌2𝐿2(Ω)+𝜌𝑡2𝐿2(Ω)𝑑𝑠.(3.28) Applying Gronwall's inequality, we obtain 𝜃(𝑡)2𝐿2(Ω)𝐶1𝜃(0)2𝐿2(Ω)+𝐶2𝑡0𝜌2𝐿2(Ω)+𝜌𝑡2𝐿2(Ω)𝑑𝑠.(3.29) Since 𝜃(0)𝐿2(Ω)𝑢(0)𝑢(0)𝐿2(Ω)+̃𝑢(0)𝑢(0)𝐿2(Ω)𝑢𝐶̃𝛾0𝐻𝛾(Ω),(3.30) we obtain the desired inequality 𝜃(𝑡)𝐿2(Ω),𝐶̃𝛾(3.31) where ̃𝛾=𝜇 if 𝜇3/2 and ̃𝛾=𝜇𝜖, 0<𝜖<1/2, if 𝜇=3/2.

Combining Lemmas 3.1 and 3.2, we obtain the following error estimates.

Theorem 3.3. Let 𝑢 and 𝑢 be the solutions of (3.1)–(3.3) and (1.2)–(1.4), respectively. Then there is a constant 𝐶(𝑢) such that 𝑢(𝑡)𝑢(𝑡)𝐿2(Ω)𝐶(𝑢)𝜇3,𝜇2,𝑢(𝑡)𝑢(𝑡)𝐿2(Ω)𝐶(𝑢)𝜇𝜖3,𝜇=21,0<𝜖<2.(3.32)

4. The Fully Discrete Variational Form

In this section, we consider a fully discrete variational formulation of (1.2). Existence and uniqueness of numerical solutions for the fully discrete variational formulation are discussed. The corresponding error estimates are also analyzed.

For the temporal discretization let 𝑘=𝑇/𝑀 for a positive integer 𝑀 and 𝑡𝑚=𝑚𝑘. Let 𝑢𝑚 be the solution of the backward Euler method defined by 𝑢𝑚+1𝑢𝑚𝑘=𝜅𝜇𝜇𝑢𝑚+1𝑢+𝑓𝑚+1(4.1) with an initial condition 𝑢0(𝑥)=𝑢0,𝑥Ω=(𝑎,𝑏)(4.2) and boundary conditions 𝑢𝑚+1(𝑎)=𝑢𝑚+1(𝑏)=0,𝑚=0,1,,𝑀1.(4.3) Then we get the fully discrete variational formulation of (1.2) to find 𝑢𝑚+1𝐻0𝜇/2(Ω) such that for all 𝑣𝐻0𝜇/2(Ω)𝑢𝑚+1𝜅,𝑣𝑘𝜇𝜇𝑢𝑚+1=𝑢,𝑣𝑘𝑓𝑚+1,𝑣+(𝑢𝑚,𝑣).(4.4) Thus a finite Galerkin solution 𝑢𝑚+1𝑉𝐻0𝜇/2(Ω) is a solution of the equation 𝑢𝑚+1,𝑣𝑘𝜅𝜇𝜇𝑢𝑚+1,𝑣𝑓𝑢=𝑘𝑚+1,𝑣+𝑢𝑚,𝑣,𝑣𝑉(4.5) with an initial condition 𝑢0=𝑢0(4.6) and boundary conditions 𝑢𝑚+1(𝑎)=𝑢𝑚+1(𝑏)=0,𝑚=0,1,,𝑀1.(4.7)

Now we prove the existence and uniqueness of solutions for (4.5) using the Brouwer fixed-point theorem.

Theorem 4.1. There exists a unique solution 𝑢𝑚+1𝑉𝐻0𝜇/2(Ω) of (4.5)–(4.7).

Proof. Let 𝐺𝑢𝑚+1=𝑢𝑚+1𝑘𝜅𝜇𝜇𝑢𝑚+1𝑢𝑘𝑓𝑚+1𝑢𝑚.(4.8) Then 𝐺(𝑣) is obviously a continuous function from 𝑉 to 𝑉. In order to show the existence of solution for 𝐺(𝑣)=0, we adopt the mathematical induction. Assume that 𝑢0,𝑢1,,𝑢𝑚 exist for 𝑚<𝑀. It follows from (1.8), Lemma 2.4, and Young's inequality that 𝑢(𝐺(𝑣),𝑣)=(𝑣,𝑣)𝑚𝜅,𝑣𝑘𝜇𝜇𝑣,𝑣𝑘(𝑓(𝑣),𝑣)𝑣2𝐿2(Ω)𝑢𝑚𝐿2(Ω)𝑣𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2||||𝑣|2𝐻0𝜇/2(Ω)𝐶𝑓𝑘𝑣2𝐿2(Ω)𝑣2𝐿2(Ω)𝑢𝑚𝐿2(Ω)𝑣𝐿2(Ω)𝐶𝑓𝑘𝑣2𝐿2(Ω)𝑣2𝐿2(Ω)12𝑢𝑚2𝐿2(Ω)+𝑣2𝐿2(Ω)𝐶𝑓𝑘𝑣2𝐿2(Ω)=12𝐶𝑓𝑘𝑣2𝐿2(Ω)12𝑢𝑚2𝐿2(Ω).(4.9) If we take sufficiently small 𝑘 so that 𝑘<1/2𝐶𝑓 and 𝑣𝐿2(Ω)>𝑢𝑚𝐿2(Ω)/(12𝐶𝑓𝑘), then the Brouwer's fixed-point theorem implies the existence of a solution.
For the proof of the uniqueness of solutions, we assume that 𝑢 and 𝑣 are two solutions of (4.5). Then we obtain (𝑢𝑣,𝜓)=𝑘𝜅𝜇(𝜇(𝑢𝑣),𝜓)+𝑘(𝑓(𝑢)𝑓(𝑣),𝜓),𝜓𝑉𝐻0𝜇/2(Ω).(4.10) Replacing 𝜓=𝑢𝑣 in the above equation and applying Lemma 2.4, we obtain 𝑢𝑣2𝐿2(Ω)𝑘𝜅𝜇|||𝜇cos𝜋2||||𝑢𝑣|𝐻0𝜇/2(Ω)+𝑘𝑓(𝑢)𝑓(𝑣)𝐿2(Ω)𝑢𝑣𝐿2(Ω)𝑘𝑓(𝑢)𝑓(𝑣)𝐿2(Ω)𝑢𝑣𝐿2(Ω)𝑘𝐶𝑙𝑢𝑣2𝐿2(Ω).(4.11) This implies 𝑢𝑣=0 since 𝑢(0)=𝑣(0).

The following theorem presents the unconditional stability for (4.4).

Theorem 4.2. The fully discrete scheme (4.4) is unconditionally stable. In fact, for any 𝑚𝑢𝑚+1𝐿2(Ω)𝑢𝐶0𝐿2(Ω).(4.12)

Proof. It follows from (1.8), Lemma 2.4, and Young's inequality that by taking 𝑣=𝑢𝑚+1 in (4.4), we obtain 𝑢0=𝑚+1,𝑢𝑚+1𝜅𝑘𝜇𝜇𝑢𝑚+1,𝑢𝑚+1𝑓𝑢𝑘𝑚+1,𝑢𝑚+1𝑢𝑚,𝑢𝑚+1𝑢𝑚+12𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2|||||𝑢𝑚+1||2𝐻0𝜇/2(Ω)𝐶𝑓𝑘𝑢𝑚+12𝐿2(Ω)𝑢𝑚𝐿2(Ω)𝑢𝑚+1𝐿2(Ω)12𝑢𝑚+12𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2|||||𝑢𝑚+1||2𝐻0𝜇/2(Ω)𝐶𝑓𝑘𝑢𝑚+12𝐿2(Ω)12𝑢𝑚2𝐿2(Ω).(4.13) Then 12𝑢𝑚+12𝐿2(Ω)12𝑢𝑚+12𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2|||||𝑢𝑚+1||2𝐻0𝜇/2(Ω)𝐶𝑓𝑘𝑢𝑚+12𝐿2(Ω)+12𝑢𝑚2𝐿2(Ω).(4.14)
Adding the above inequality from 𝑚=0 to 𝑚, we obtain 12𝐶𝑓𝑘𝑢𝑚+12𝐿2(Ω)𝑢02𝐿2(Ω)+2𝐶𝑓𝑘𝑚𝑗=1𝑢𝑗2𝐿2(Ω).(4.15) Applying the discrete Gronwall's inequality with sufficiently small 𝑘 such that 𝑘<1/2𝐶𝑓, we obtain the desired result.

The following theorem is an error estimate for the fully discrete problem (4.4).

Theorem 4.3. Let 𝑢 be the exact solution of (1.2) and let 𝑢𝑚 be the solution of (4.4). Then there is a constant 𝐶 such that 𝑢𝑡𝑚𝑢𝑚𝐿2(Ω)𝐶𝑘.(4.16)

Proof. Let 𝑒𝑚=𝑢(𝑡𝑚)𝑢𝑚 be the error at 𝑡𝑚. It follows from (1.2) and (4.4) that for any 𝑣𝐻0𝜇/2(Ω)𝑒𝑚+1𝜅,𝑣𝑘𝜇𝜇𝑒𝑚+1𝑓𝑢𝑡,𝑣=𝑘𝑚+1𝑢𝑓𝑚+1,𝑣+(𝑒𝑚,𝑣)+𝑘𝑟𝑚+1,,𝑣(4.17) where 𝑟=𝑂(𝑘). Taking 𝑣=𝑒𝑚+1, 𝑒𝑚+12𝐿2(Ω)𝑒𝑚+12𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2|||||𝑒𝑚+1||2𝐻0𝜇/2(Ω)𝑘𝑓(𝑢(𝑡𝑚+1))𝑓(𝑢𝑚+1)𝐿2(Ω)𝑒𝑚+1𝐿2(Ω)+𝑒𝑚𝐿2(Ω)𝑒𝑚+1𝐿2(Ω)+𝑘𝑟𝑚+1𝐿2(Ω)𝑒𝑚+1𝐿2(Ω).(4.18) Applying the locally Lipschitz continuity of 𝑓 and Young's inequality, we obtain 𝑒𝑚+12𝐿2(Ω)𝑘𝐶𝑙𝑒𝑚+12𝐿2(Ω)+𝑒𝑚𝐿2(Ω)𝑒𝑚+1𝐿2(Ω)+𝑘𝑟𝑚+1𝐿2(Ω)𝑒𝑚+1𝐿2(Ω)𝑘𝐶𝑙𝑒𝑚+12𝐿2(Ω)+𝜀1𝑒𝑚2𝐿2(Ω)+14𝜀1𝑒𝑚+12𝐿2(Ω)+𝜀2𝑘𝑟𝑚+12𝐿2(Ω)+14𝜀2𝑒𝑚+12𝐿2(Ω).(4.19) That is, 114𝜀114𝜀2𝑒𝑚+12𝐿2(Ω)𝑘𝐶𝑙𝑒𝑚+12𝐿2(Ω)+𝜀1𝑒𝑚2𝐿2(Ω)+𝜀2𝑘𝑟𝑚+12𝐿2(Ω).(4.20) Denoting 𝜀0=11/4𝜀11/4𝜀2 and adding the above equation from 𝑚=0 to 𝑚, we obtain 𝜀0𝑘𝐶𝑙𝑒𝑚+12𝐿2(Ω)𝜀1𝑒02𝐿2(Ω)+𝑘𝐶𝑙+𝜀1𝜀0𝑚𝑖=1𝑒𝑖2𝐿2(Ω)+𝜀2𝑚+1𝑖=1𝑘𝑟𝑖2𝐿2(Ω).(4.21) Applying the discrete Gronwall's inequality with sufficiently small 𝑘 such that (𝜀0𝜀1)/𝐶𝑙<𝑘<𝜀0/𝐶𝑙, we obtain the desired result since 𝑚+1𝑖=1𝑘𝑟𝑖𝐿2(Ω)𝐶𝑘 and 𝑒0𝐿2(Ω)=𝑢(0)𝑢0𝐿2(Ω)=0.

As in the previous section, denote 𝜃𝑚+1=𝑢𝑚+1̃𝑢𝑚+1 and 𝜌𝑚+1=̃𝑢𝑚+1𝑢(𝑡𝑚+1). Here ̃𝑢𝑚+1 is the elliptic projection of 𝑢(𝑡𝑚+1) defined in (3.13). Then 𝑒𝑚+1=𝜃𝑚+1+𝜌𝑚+1.(4.22)

Theorem 4.4. Let 𝑢 be the exact solution of (1.2)–(1.4) and let {𝑢𝑚}𝑀𝑚=0 be the solution of (4.5)–(4.7). Then when 𝜇3/2𝑢(𝑡𝑚+1)𝑢𝑚+1𝐿2(Ω)𝐶𝑘+𝐶𝜇𝑢𝑡𝑚+1𝐻𝜇(Ω)(4.23) and when 𝜇=3/2,0<𝜖<1/2, 𝑢(𝑡𝑚+1)𝑢𝑚+1𝐿2(Ω)𝐶𝑘+𝐶𝜇𝜖𝑢(𝑡𝑚+1)𝐻𝜇𝜖(Ω).(4.24)

Proof. Since we know the estimates on 𝜌 from Lemma 3.1, we have only to show boundedness of 𝜃𝑚+1. Using the property (3.13), we obtain for 𝑣𝑉𝜃𝑚+1𝜅,𝑣𝑘𝜇𝜇𝜃𝑚+1𝑓𝑢,𝑣=𝑘𝑚+1𝑢𝑡𝑓𝑚+1+𝑢,𝑣𝑚𝑡𝑢𝑚,𝑣𝑘𝑟𝑚+1𝜌,𝑣𝑚+1,,𝑣(4.25) where 𝑟=𝑂(𝑘).
Taking 𝑣=𝜃𝑚+1 and applying Lemma 2.4, the locally Lipschitz continuity of 𝑓, Young's inequality, and the triangle inequality, we obtain 𝜃𝑚+12𝐿2(Ω)𝜃𝑚+12𝐿2(Ω)+𝑘𝜅𝜇|||𝜇cos𝜋2|||||𝜃𝑚+1||2𝐻0𝜇/2(Ω)𝑘𝑓(𝑢𝑚+1)𝑓(𝑢(𝑡𝑚+1))𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝑒𝑚𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝑘𝑟𝑚+1𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝜌𝑚+1𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)𝑘𝐶𝑙𝑒𝑚+1𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝜃𝑚𝐿2(Ω)+𝜌𝑚𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝑘𝑟𝑚+1𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)+𝜌𝑚+1𝐿2(Ω)𝜃𝑚+1𝐿2(Ω)𝑘𝐶𝑙11+4𝜀6𝜃𝑚+12𝐿2(Ω)+14𝜀3+14𝜀4+14𝜀5+14𝜀6𝜃𝑚+12𝐿2(Ω)+𝜀3𝜃𝑚2𝐿2(Ω)+𝜀4𝑘𝑟𝑚+12𝐿2(Ω)+𝜀5𝜌𝑚2𝐿2(Ω)+1+𝑘𝐶𝑙𝜀6𝜌𝑚+12𝐿2(Ω).(4.26) This implies that 114𝜀314𝜀414𝜀514𝜀6𝜃𝑚+12𝐿2(Ω)𝑘𝐶𝑙11+4𝜀6𝜃𝑚+12𝐿2(Ω)+𝜀3𝜃𝑚2𝐿2(Ω)+𝜀4𝑘𝑟𝑚+12𝐿2(Ω)+𝜀5𝜌𝑚2𝐿2(Ω)+1+𝑘𝐶𝑙𝜀6𝜌𝑚+12𝐿2(Ω).(4.27) Denote 𝜀7=11/4𝜀31/4𝜀41/4𝜀51/4𝜀6 and 𝜀8=1+1/4𝜀6. Then adding the above inequality from 𝑚=0 to 𝑚, we obtain 𝜀7𝑘𝐶𝑙𝜀8𝜃𝑚+12𝐿2(Ω)𝜀3𝜃02𝐿2(Ω)+𝑘𝐶𝑙𝜀8+𝜀3𝜀7𝑚𝑖=1𝜃𝑖2𝐿2(Ω)+𝜀4𝑚+1𝑖=1𝑘𝑟𝑖2𝐿2(Ω)+𝜀5𝑚𝑖=0𝜌𝑖2𝐿2(Ω)+1+𝑘𝐶𝑙𝜀6𝑚+1𝑖=1𝜌𝑖2𝐿2(Ω).(4.28) Applying the discrete Gronwall's inequality with sufficiently small 𝑘 such that (𝜀7𝜀3)/𝜀8𝐶𝑙<𝑘<𝜀7/𝐶𝑙𝜀8, 𝜃𝑚+12𝐿2(Ω)𝐶1𝜃02𝐿2(Ω)+𝐶2𝑚+1𝑖=1𝑘𝑟𝑖2𝐿2(Ω)+𝐶3𝑚+1𝑖=0𝜌𝑖2𝐿2(Ω).(4.29) Also, using Lemma 3.1 and the initial conditions (1.3) and (4.6), we obtain 𝜃0𝐿2(Ω)𝑢0𝑢(0)𝐿2(Ω)+̃𝑢0𝑢(0)𝐿2(Ω)𝑢𝐶̃𝛾0𝐻𝛾(Ω).(4.30) Since 𝑚+1𝑖=1𝑘𝑟𝑖𝐿2(Ω)𝐶𝑘, we get 𝜃𝑚+1𝐿2(Ω),𝐶𝑘+𝐶(𝑢)̃𝛾(4.31) where ̃𝛾=𝜇 if 𝜇3/2 and ̃𝛾=𝜇𝜖, 0<𝜖<1/2, if 𝜇=3/2. Thus we obtain the desired result.

5. Numerical Experiments

In this section, we present numerical results for the Galerkin approximations which supports the theoretical analysis discussed in the previous section.

Let 𝑆 denote a uniform partition of Ω and let 𝑉 denote the space of continuous piecewise linear functions defined on 𝑆. In order to implement the Galerkin finite element approximation, we adapt finite element discretization on the spatial axis and the backward Euler finite difference scheme along the temporal axis. We associate shape functions of space 𝑉 with the standard basis of the functions on the uniform interval with length .

Example 5.1. We first consider a space fractional linear diffusion equation: 𝜕𝑢(𝑥,𝑡)𝜕𝑡=𝜇𝑢(𝑥,𝑡)+2𝑡𝑡2𝑡+1𝑢(𝑥,𝑡)2×𝑥+12𝜇+(1𝑥)2𝜇Γ6𝑥(3𝜇)3𝜇+(1𝑥)3𝜇Γ+𝑥(4𝜇)124𝜇+(1𝑥)4𝜇Γ(5𝜇)(5.1) with an initial condition 𝑢(𝑥,0)=𝑥2(1𝑥)2[],𝑥0,1(5.2) and boundary conditions 𝑢(0,𝑡)=𝑢(1,𝑡)=0.(5.3) In this case, the exact solution is 𝑡𝑢(𝑥,𝑡)=2𝑥+12(1𝑥)2.(5.4)

Tables 1, 2, and 4 show the order of convergence and 𝐿2-error between the exact solution and the Galerkin approximate solution of the fully discrete backward Euler method for (5.1) when 𝜇=1.6, 𝜇=1.8 and 𝜇=1.5, respectively. For numerical computation, the temporal step size 𝑘=0.001 is used in all three cases. Table 3 shows 𝐿2-errors and orders of convergence for the Galerkin approximate solution when 𝜇=1.8 and the spatial step size =0.0625.

According to Tables 13, we may find the order of convergence of 𝑂(𝑘+𝜇) for this linear fractional diffusion problem (5.1)–(5.3) when 𝜇3/2. Furthermore, Table 4 shows orders of numerical convergence for the problem when 𝜇=3/2, where we may see that the order of convergence is of 𝑂(𝑘+𝜇𝜖), 0<𝜖<1/2. It follows from Tables 14 that numerical computations confirm the theoretical results.

We plot the exact solution and approximate solutions obtained by the backward Euler Galerkin method using =1/32 and 𝑘=1/1000 for (5.1) with 𝜇=1.6 and 𝜇=1.8. Figure 1 shows the contour plots of an exact solution and numerical solutions at 𝑡=1, and Figure 2 shows log-log graph for the order of convergence.

Example 5.2. We consider a space fractional diffusion equation with a nonlinear Fisher type source term which is described as 𝜕𝑢(𝑥,𝑡)𝜕𝑡=𝜅𝜇𝜇𝑢(𝑥,𝑡)+𝜆𝑢(𝑥,𝑡)(1𝛽𝑢(𝑥,𝑡))(5.5) with an initial condition 𝑢(𝑥,0)=𝑢0(𝑥)(5.6) and boundary conditions 𝑢(1,𝑡)=𝑢(1,𝑡)=0.(5.7) In fact, we will consider the case of 𝜅𝜇=0.1, 𝛽=1 in (5.5) with an initial condition 𝑢0𝑒(𝑥)=10𝑥𝑒,𝑥0,10𝑥,𝑥<0.(5.8)

For numerical computations, we have to take care of the nonlinear term 𝑓(𝑢)=𝜆𝑢(1𝛽𝑢). This gives a complicated nonlinear matrix. In order to avoid the difficulty of solving nonlinear system, we adopted a linearized method replacing 𝜆𝑢𝑛+1(1𝛽𝑢𝑛+1) by 𝜆𝑢𝑛+1(1𝛽𝑢𝑛). Figure 3 shows contour plots of numerical solutions at 𝑡=1 for (5.5)–(5.8) with 𝜆=0.25. For numerical computations, step sizes =0.01 and 𝑘=0.005 are used. From the numerical results we may find that numerical solutions converge to the solution of classical diffusion equation as 𝜇 approaches to 2.

Example 5.3. We now consider (5.5) with 𝜅𝜇=0.1, 𝛽=1 and boundary conditions lim|𝑥|𝑢(𝑥,𝑡)=0.(5.9) We will consider an initial condition with a sharp peak in the middle as 𝑢0(𝑥)=sech2(10𝑥)(5.10) and an initial condition with a flat roof in the middle as 𝑢0𝑒(𝑥)=10(𝑥1)𝑒,𝑥>1,1,1<𝑥1,10(𝑥+1),𝑥1.(5.11)

Tang and Weber [30] have obtained computational solutions for (5.5) with initial conditions (5.10) and (5.11) using a Petrov-Galerkin method when (5.5) is a classical diffusion problem. We obtain computational results using the method as in Example 5.2. Figure 4 shows contour plots of numerical solutions at 𝑡=1 for (5.5) with an initial condition (5.10) when Ω=(2,2) and 𝜆=0.25. Figure 5 shows also contour plots of numerical solutions at 𝑡=4 for (5.5) and (5.10) when Ω=(4,4) and 𝜆=1. In both cases, step sizes =0.01 and 𝑘=0.005 are used for computation. According to Figures 4 and 5, we may see that the diffusivity depends on 𝜇 but it is far less than that of the classical solution. That is, the fractional diffusion problem keeps the peak in the middle for longer time than the classical one does.

Figure 6 shows contour plots of numerical solutions for (5.5) with an initial condition (5.10) when 𝜇=1.8, Ω=(2,2) and 𝜆=1. In this case, step sizes =0.01 and 𝑘=0.005 are also used for computation. But the period of time is from 𝑡=0 to 𝑡=5. According to Figure 6, we may see that the peak goes down rapidly for a short time, and it begins to go up after the contour arrives at the lowest level.

Figure 7 shows contour plots of numerical solutions at 𝑡=1 for (5.5) with an initial condition (5.11) when Ω=(4,4) and 𝜆=0.25. In this case, step sizes =0.01 and 𝑘=0.005 are also used for computation. According to Figure 7, we may find that the fractional diffusion problem keeps the flat roof in the middle for longer time than the classical one does.

6. Concluding Remarks

Galerkin finite element methods are considered for the space fractional diffusion equation with a nonlinear source term. We have derived the variational formula of the semidiscrete scheme by using the Galerkin finite element method in space. We showed existence and stability of solutions for the semidiscrete scheme. Furthermore, we derived the fully time-space discrete variational formulation using the backward Euler method. Existence and uniqueness of solutions for the fully discrete Galerkin method have been discussed. Also we proved that the scheme is unconditionally stable, and it has the order of convergence of )𝑂(𝑘+̃𝛾, where ̃𝛾 is a constant depending on the order of fractional derivative. Numerical computations confirm the theoretical results discussed in the previous section for the problem with a linear source term. For the fractional diffusion problem with a nonlinear source term, we may find that the diffusivity depends on the order of fractional derivative, and numerical solutions of fractional order problems are less diffusive than the solution of a classical diffusion problem.

Acknowledgment

The authors would like to express sincere thanks to the referee for their invaluable comments.