A Fully Discrete Galerkin Method for a Nonlinear Space-Fractional Diffusion Equation

Abstract

The spatial transport process in fractal media is generally anomalous. The space-fractional advection-diffusion equation can be used to characterize such a process. In this paper, a fully discrete scheme is given for a type of nonlinear space-fractional anomalous advection-diffusion equation. In the spatial direction, we use the finite element method, and in the temporal direction, we use the modified Crank-Nicolson approximation. Here the fractional derivative indicates the Caputo derivative. The error estimate for the fully discrete scheme is derived. And the numerical examples are also included which are in line with the theoretical analysis.

1. Introduction

The normal diffusive motion is modeled to describe the standard Brownian motion. The relation between the flow and the divergence of the particle displacement represents 𝐽(𝑥,𝑡)=−𝑎𝜕𝑐𝜕𝑥+𝑏𝑐,(1.1) where 𝐽 is the diffusive flow. Inserting the above equation into the equation of mass conservation 𝜕𝐽𝜕𝑥=−𝜕𝜕𝑐𝑡,(1.2) we obtain the standard convection-diffusion equation. From the viewpoint of physics, it means that during the method of time random walkers, the overall particle displacement up to time 𝑡 can be represented as a sum of independent random steps, in the case that both the mean-squared displacement per step and the mean time needed to perform a step are finite. The measured variance growth in the direction of flow of tracer plumes is typically at a Fickian rate, ⟨(𝑐−𝑐)²⟩∼𝑡.

The transport process in fractal media cannot be described with the normal diffusion. The process is nonlocal and it does not follow the classical Fickian law. It depicts a particle in spreading tracer cloud which has a standard deviation, and which grows like 𝑡^2𝛼 for some 0<𝛼<1, excluding the Fickian case 𝛼=1/2. The description of anomalous diffusion means that the measure variance growth in the direction of flow has a deviation from the Fickian case, it follows the super-Fickian rate ⟨(𝑐−𝑐)²⟩∼𝑡^2𝛼 when 𝛼>1/2, or does the subdiffusion rate ⟨(𝑐−𝑐)²⟩∼𝑡^2𝛼 if 0<𝛼<1/2. With the help of the continuous time random walk and the Fourier transform, the governing equation with space fractional derivative can be derived as follows 𝜕𝑢𝜕𝑡=𝐷𝑎(𝑢)^𝑎𝐷^𝛽𝑥𝑢+𝑏(𝑢)𝐷𝑢+𝑓(𝑥,𝑡,𝑢),0<𝛽<1,(1.3) where 𝐷 denotes integer derivative respect to 𝑥, and 𝐷^𝛽 is fractional derivative. There are some authors studying the spacial anomalous diffusion equation in theoretical analysis and numerical simulations [1–10]. Now the fractional anomalous diffusion becomes a hot topic because of its widely applications in the evolution of various dynamical systems under the influence of stochastic forces. For example, it is a well-suited tool for the description of anomalous transport processes in both absence and presence of external velocities or force fields. Since the groundwater velocities span many orders of magnitude and give rise to diffusion-like dispersion (a term that combines molecular diffusion and hydrodynamic dispersion), the fractional diffusion is an important process in hydrogeology. It can be used to describe the systems with reactions and diffusions across a wide range of applications including nerve cell signaling, animal coat patterns, population dispersal, and chemical waves. In general, fractional anomalous diffusions have numerous applications in statistical physics, biophysics, chemistry, hydrogeology, and biology [4, 11–20].

In this paper, we mainly study one kind of typical nonlinear space-fractional partial differential equations by using the finite element method, which reads in the following form: 𝜕𝑢𝜕𝑡=𝐷𝑎(𝑢)^𝑎𝐷^𝛽𝑥𝑢𝑢+∣𝑏(_𝑡=0𝑢)𝐷𝑢+𝑓(𝑥,𝑡,𝑢),𝑥∈Ω,𝑡∈(0,𝑇],𝑢∣_𝜕Ω=𝜑(𝑥),𝑥∈Ω,=𝑔,𝑡∈(0,𝑇],(1.4) where Ω is a spacial domain with boundary 𝜕Ω, 𝐷^𝛽 is the 𝛽th (0<𝛽<1) order fractional derivative with respect to the space variable 𝑥 in the Caputo sense (which will be introduced later on),𝑎,𝑏,𝑓 are functions of𝑥,𝑡,𝑢,𝜑 and 𝑔 are known functions which satisfy the conditions requested by the theorem of error estimations.

The rest of this paper is constructed as follows. In Section 2 the fractional integral, fractional derivative, and the fractional derivative spaces are introduced. The error estimates of the finite element approximation for (1.4) are studied in Section 3, and in Section 4, numerical examples are taken to verify the theoretical results derived in Section 3.

2. Fractional Derivative Space

In this section, we firstly introduce the fractional integral (or Riemann-Liouville integral), the Caputo fractional derivative, and their corresponding fractional derivative space.

Definition 2.1. The 𝛼th order left and right Riemann-Liouville integrals of function 𝑢(𝑥) are defined as follows _𝑎𝐼^𝛼𝑥𝑢(𝑥)=1Γ(𝛼)^𝑥𝑎(𝑥−𝑠)^𝛼−1𝑢(𝑠)𝑑𝑠,_𝑥𝐼^𝛼𝑏𝑢(𝑥)=1Γ(𝛼)^𝑏𝑥(𝑠−𝑥)^𝛼−1𝑢(𝑠)𝑑𝑠,(2.1) where 𝛼>0, and Γ(⋅) is the Gamma function.

Definition 2.2. The 𝛼th order Caputo derivative of function 𝑢(𝑥) is defined as, _𝑎𝐷^𝛼𝑥𝑢(𝑥)=_𝑎𝐼^{𝑥𝑛−𝛼}𝑑^𝑛𝑢(𝑥)𝑑𝑥^𝑛,𝑛−1<𝛼<𝑛∈𝑍⁺,_𝑥𝐷^𝛼𝑏𝑢(𝑥)=(−1)^𝑛𝑥𝐼^{𝑏𝑛−𝛼}𝑑^𝑛𝑢(𝑥)𝑑𝑥^𝑛,𝑛−1<𝛼<𝑛∈𝑍⁺.(2.2) The 𝛼th order Riemann-Liouville derivative of function 𝑢(𝑥) is defined by changing the order of integration and differentiation.

Lemma 2.3 (see [8]). If 𝑢(0)=𝑢′(0)=⋯=𝑢^(𝑛−1)(0)=0, then the Caputo fractional derivative is equal to the Riemann-Liouville derivative.

Definition 2.4. The fractional derivative space 𝐽^𝛼(Ω) is defined as follows: 𝐽^𝛼(Ω)=𝑢∈𝐿²(Ω)∶_𝑎𝐷^𝛼𝑥𝑢∈𝐿²(Ω),𝑛−1≤𝛼<𝑛,(2.3) endowed with the seminorm |𝑢|_𝐽^𝛼=‖‖_𝑎𝐷^𝛼𝑥𝑢‖‖_𝐿2(Ω),(2.4) and the norm ‖𝑢‖_𝐽^𝛼=|𝑢|^2𝐽_𝛼+_{[𝛼]𝑘≤}‖‖𝐷^𝑘𝑢‖‖²^1/2.(2.5)

Let 𝐽^𝛼0(Ω) denote the closure of 𝐶^∞0(Ω) with respect to the above norm and seminorm.

Definition 2.5. Define the seminorm |𝑢|_𝐻^𝛼=‖|𝑖𝑤|^𝛼𝐹(𝑢)‖_𝐿²_(Ω),(2.6) and the norm ‖𝑢‖_𝐻^𝛼=|𝑢|^2𝐻_𝛼+_{[𝛼]𝑘≤}‖‖𝐷^𝑘𝑢‖‖²^1/2,(2.7) where 𝑖 is the imaginary unit, and 𝐹 is the Fourier transform, and which can define another fractional derivative space 𝐻^𝛼(Ω).

Let 𝐻^𝛼0(Ω) denote the closure of 𝐶^∞0(Ω) with respect to the norm and seminorm.

Definition 2.6. The fractional space 𝐽^𝛼𝑠(Ω) is defined below 𝐽^𝛼𝑠(Ω)=𝑢∈𝐿²(Ω)∶_𝑎𝐷^𝛼𝑥𝑢∈𝐿²(Ω),_𝑥𝐷^𝛼𝑏𝑢∈𝐿²(Ω),𝑛−1≤𝛼<𝑛,(2.8) endowed with the seminorm |𝑢|_𝐽^𝛼𝑠=|||_𝑎𝐷^𝛼𝑥𝑢,_𝑥𝐷^𝛼𝑏𝑢^1/2|||_𝐿2(Ω),(2.9) and the norm ‖𝑢‖_𝐽^𝛼𝑠=_{[𝛼]𝑘≤}‖‖𝐷^𝑘𝑢‖‖²+|𝑢|^2𝐽_𝛼𝑠^1/2.(2.10)

Theorem 2.7 (see [3, 6]). 𝐽^𝛼𝑠, 𝐽^𝛼, and 𝐻^𝛼 are equal with equivalent seminorm and norm.
The following are some useful results.

Lemma 2.8 (see [3]). For 𝑢∈𝐽^𝛼0(Ω), 0<𝛽<𝛼, then _𝑎𝐷^𝛼𝑥𝑢(𝑥)=_𝑎𝐷^{𝑥𝑎𝛼−𝛽}𝐷^𝛽𝑥𝑢.(2.11)

Lemma 2.9 (see [2]). For 𝑢∈𝐻^𝛼0(Ω), one has ‖𝑢‖_𝐿²_(Ω)≤𝑐|𝑢|_𝐻^𝛼0.(2.12) For 0<𝛽<𝛼, |𝑢|_𝐻^𝛽0_(Ω)≤𝑐|𝑢|_𝐻^𝛼0.(2.13)

Since 𝐽^𝛼𝑠, 𝐽^𝛼, and 𝐻^𝛼 are equal with equivalent seminorm and norm, the norms with each space which will be used following are without distinction, and the notations are used seminorm |⋅|_𝛼 and norm ‖⋅‖_𝛼.

3. Finite Element Approximation

Let Ω=[𝑎,𝑏], and 0≤𝛽<1. Define 𝛼=(1+𝛽)/2. In this section, we will formulate a fully discrete Galerkin finite element method for a type of nonlinear anomalous diffusion equation as follows.

Problem 1 (Nonlinear spacial anomalous diffusion equation). We consider equations of the form 𝜕𝑢𝜕𝑡=𝐷𝑎(𝑢)^𝑎𝐷𝑢(^𝛽𝑥𝑥𝑢,𝑡)𝑢+=(𝑏𝑥𝜙,(𝑢0())𝑥𝐷=,𝑡𝑢)𝑔,+(𝑥𝑓𝑥),(∈𝑥,𝑥𝜕𝑡,Ω∈𝑢)×,(0,(𝑥𝑇,]𝑡,)∈Ω×(0,𝑇],Ω.(3.1) We always assume that 0<𝑚<𝑎(𝑢)<𝑀,0<𝑚<𝑏(𝑢)<𝑀,0<𝑚<𝑓(𝑢)<𝑀.(3.2)

The algorithm and analysis in this paper are applicable for a large class of linear and nonlinear functions (including polynomials and exponentials) in the unknown variables. Throughout the paper, we assume the following mild Lipschitz continuity conditions on  𝑎,𝑏,  and  𝑓: there exist positive constants 𝐿 and 𝑐 such that for 𝑥∈Ω,𝑡∈(0,𝑇], and 𝑠,𝑡∈𝑅, ||𝑎(𝑥,𝑡,𝑠)−𝑎(𝑥,𝑡,𝑟)||≤𝐿|𝑠−𝑟|,(3.3)||𝑏(𝑥,𝑡,𝑠)−𝑏(𝑥,𝑡,𝑟)||≤𝐿|𝑠−𝑟|,(3.4)||𝑓(𝑥,𝑡,𝑠)−𝑓(𝑥,𝑡,𝑟)||≤𝐿|𝑠−𝑟|.(3.5)

In order to derive a variational form of Problem 1, we suppose that 𝑢 is a sufficiently smooth solution of Problem 1. Multiplying an arbitrary 𝑣∈𝐻^𝛼0(Ω) in both sides yields _Ω𝜕𝑢𝜕𝑡𝑣𝑑𝑥=_Ω𝐷𝑎(𝑢)_𝑎𝐷^𝛽𝑥𝑢𝑣𝑑𝑥+_Ω𝑏(𝑢)𝐷𝑢𝑣𝑑𝑥+_Ω𝑓(𝑥,𝑡,𝑢)𝑣𝑑𝑥.(3.6)

Rewriting the above expression yields _Ω𝜕𝑢𝜕𝑡𝑣𝑑𝑥+_Ω𝑎(𝑢)_𝑎𝐷^𝛽𝑥𝑢𝐷𝑣𝑑𝑥−_Ω𝑏(𝑢)𝐷𝑢𝑣𝑑𝑥=_Ω𝑓(𝑥,𝑡,𝑢)𝑣𝑑𝑥.(3.7)

We define the associated bilinear form 𝐴∶𝐽^𝛼0(Ω)×𝐽^𝛼0(Ω)→𝑅 as 𝐴(𝑢,𝑣)=𝑎(𝑢)_𝑎𝐷^𝛽𝑥𝑢,𝐷𝑣−(𝑏(𝑢)𝐷𝑢,𝑣),(3.8) where (⋅,⋅) denotes the inner product on 𝐿²(Ω) and 𝐽^𝛼0(Ω).

For given 𝑓∈𝐽^−𝛼(Ω), we define the associated function 𝐹∶𝐽^𝛼0(Ω)→𝑅 as 𝐹(𝑣)=⟨𝑓,𝑣⟩.(3.9)

Definition 3.1. A function 𝑢∈𝐽^𝛼0(Ω) is a variational solution of Problem 1 provided that 𝜕𝜕𝑢𝑡,𝑣+𝐴(𝑢,𝑣)=𝐹(𝑣),∀𝑣∈𝐽^𝛼0(Ω).(3.10)

Now we are ready to describe a fully discrete Galerkin finite element method to solve nonlinear Problem 1. In our new scheme, the finite element trial and test spaces for Problem 1 are chosen to be same.

For a positive integer 𝑁, let ∏^𝑡={𝑡_𝑛}^𝑁𝑛=0 be a uniform partition of the time interval (0,𝑇] such that 𝑡_𝑛=𝑛𝜏, where 𝜏=𝑇/𝑁, and let 𝑡_𝑛−1/2=𝑡_𝑛−𝜏/2. Throughout the paper, we use the following notation for a function 𝜙: 𝜙^𝑛=𝜙𝑡_𝑛,𝜕_𝑡𝜙^𝑛=𝜙^𝑛−𝜙^𝑛−1𝜏,𝜙^𝑛=𝜙^𝑛+𝜙^𝑛−12,𝜙^𝑛=3𝜙^𝑛−1−𝜙^𝑛−22.(3.11)

Let 𝒦_ℎ={𝐾} be a partition of spatial domain Ω. Define ℎ_𝑘 as the diameter of the element 𝐾 and ℎ=max_{𝐾∈𝒦ℎ}ℎ_𝐾. And let 𝑆_ℎ be a finite element space 𝑆_ℎ=𝑣∈𝐻^𝛼0(Ω)∶𝑣∣_𝐾∈𝑃_𝑟−1(𝐾),𝐾∈𝒦_ℎ,(3.12) where 𝑃_𝑟−1(𝐾) is the set of polynomials of degree 𝑟−1 on a given domain 𝐾. And the functions in 𝑆_ℎ are continuous on Ω. Our fully discrete quadrature scheme to solve Problem 1 is to find 𝑢_ℎ: for 𝑣∈𝑆_ℎ such that 𝜕_𝑡𝑢^𝑛ℎ,𝑣+𝑎𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝑢^𝑛ℎ,𝐷𝑣−𝑏𝑢^𝑛ℎ̃𝐷𝑢^𝑛ℎ,𝑣=𝑓𝑢^𝑛ℎ̃,𝑣.(3.13)

The linear systems in the above equation requires selecting the value of 𝑢^0ℎ and 𝑢^1ℎ. Given 𝑢^0ℎ depending on the initial data 𝑔(𝑥), we select 𝑢^1ℎ by solving the following predictor-corrector linear systems: 𝑢^ℎ1,0−𝑢^0ℎ𝜏,𝑣+𝑎𝑢^0ℎ_𝑎𝐷^𝛽𝑥𝑢^ℎ1,0+𝑢^0ℎ2,𝐷𝑣−𝑏𝑢^0ℎ𝐷𝑢^ℎ1,0+𝑢^0ℎ𝑢2,𝑣=𝑓𝑢^0ℎ,𝑣_1ℎ,−𝑢^0ℎ𝜏,𝑣+𝑎𝑢^ℎ1,0+𝑢^0ℎ2^𝑎𝐷^𝛽𝑥𝑢^1ℎ+𝑢^0ℎ2,𝐷𝑣−𝑏𝑢^ℎ1,0+𝑢^0ℎ2𝐷𝑢^1ℎ+𝑢^0ℎ=𝑓𝑢2,𝑣_ℎ1,0+𝑢^0ℎ2,𝑣.(3.14)

Lemma 3.2. For 𝑢,𝑣,𝑤∈𝐽^𝛼𝑠,0(Ω),0<𝑚≤𝑎(𝑢)≤𝑀,𝛼=(1+𝛽)/2, there exist constants 𝛾₁,𝛾₂ such that 𝑎(𝑢)^𝑎𝐷^𝛽𝑥𝑢,𝐷𝑣≤𝛾₁‖𝑢‖_𝛼⋅‖𝑣‖_𝛼,𝑎(𝑤)_𝑎𝐷^𝛽𝑥𝑣,𝐷𝑣≥𝛾₂‖𝑣‖^2𝛼.(3.15)

Proof. With the assumption of 𝑎(𝑢) in (3.3) and the property of dual space 𝑎(𝑤)^𝑎𝐷^𝛽𝑥𝑢,𝐷𝑣≤‖‖𝑎(𝑤)_𝑎𝐷^𝛽𝑥𝑢‖‖_1−𝛼⋅‖𝐷𝑣‖_{−(1−𝛼)}𝑎(𝑤)_𝑎≤𝑀𝑐‖𝑢‖_{1−𝛼+𝛽}⋅‖𝑣‖_{−(1−𝛼)+1}≤𝛾₁‖𝑢‖_𝛼⋅‖𝑣‖_𝛼,𝐷^𝛽𝑥𝑣,𝐷𝑣==−−𝐷_𝑎𝑎(𝑤)_𝑎𝐷^𝛽𝑥𝑣,𝑣𝐷^{𝑥(1−𝛽)/2}𝑎(𝑤)_𝑎𝐷^𝛽𝑥𝑣,_𝑥𝐷^{𝑏(1+𝛽)/2}𝑣≥𝑚|𝑣|^2𝐽_𝛼𝑠≥𝛾₂‖𝑣‖^2𝛼.(3.16)

Lemma 3.3 (see [2]). For Ω⊂𝑅^𝑛,𝛼>𝑛/4,𝑣,𝑤∈𝐻^𝛼0(Ω),𝜀>0, one has (𝑣𝑏(𝑤),∇𝑣)≤𝑐₀(𝑞𝜀)^−𝑝/𝑞𝑝‖∇𝑏(𝑤)‖^𝑝⋅‖𝑣‖²+𝜀‖𝑣‖^2𝛼,(3.17) where 𝑝=4𝛼/(4𝛼−𝑛),𝑞=4𝛼/𝑛.

Theorem 3.4. Let 𝑢^𝑛ℎ be bounded, then for a sufficiently small step 𝜏, there exists a unique solution 𝑢^𝑛ℎ∈𝑆_ℎ satisfying scheme (3.13).

Proof. As scheme represents a finite system of problem, the continuity and coercivity of (𝑢^𝑛ℎ,𝜔^𝑛ℎ)/𝜏+𝐴(𝑢^𝑛ℎ,𝜔^𝑛ℎ) is the sufficient and essential condition for the existence and uniqueness of 𝑢^𝑛ℎ. Let 𝑣=𝑢^𝑛ℎ,𝑤=𝜔^𝑛ℎ, then (𝑣,𝑣)𝜏+𝐴(𝑣,𝑣)=(𝑣𝜏,𝑣)+𝑎(𝑤)^𝑎≥‖𝑣‖₂𝐷^𝛽𝑥𝑣,𝐷𝑣−(𝑏(𝑤)𝐷𝑣,𝑣)𝜏+𝛾²‖𝑣‖_𝛼−𝑐₀‖𝐷𝑏(𝑤)‖²‖𝑣‖²−𝜀‖𝑣‖^2𝛼=𝛾₂−𝜀‖𝑣‖^2𝛼+𝜏⁻¹−𝑐₀‖𝐷𝑏(𝑤)‖²‖𝑣‖²≥𝑐‖𝑣‖^2𝛼.(3.18) For the chosen sufficiently small 𝜏, the above inequality holds. (𝑣,𝑤)𝜏+𝐴(𝑣,𝑤)=(𝑣,𝜏𝑤)+𝑎(𝑢)^𝑎≤‖𝑢‖⋅‖𝑤‖𝐷^𝛽𝑥𝑣,𝐷𝑤+(𝐷𝑏(𝑢)𝑣,𝐷𝑤)𝜏+𝛾¹≤‖𝑢‖⋅‖𝑤‖‖𝑣‖_𝛼‖𝑤‖_𝛼+‖𝑣‖⋅‖𝐷(𝑏(𝑢)𝑤)‖𝜏+𝛾¹‖𝑣‖_𝛼‖𝑤‖_𝛼+𝑀‖𝑣‖⋅‖𝑤‖≤𝑐‖𝑣‖_𝛼ℎ‖𝑤‖_𝛼.(3.19) Hence, the scheme (3.13) is uniquely solvable for 𝑢^𝑛ℎ.
Let 𝜌^𝑛=𝑃_ℎ𝑢^𝑛−𝑢^𝑛, and 𝜃^𝑛=𝑢^𝑛ℎ−𝑃_ℎ𝑢^𝑛, then 𝑢^𝑛ℎ−𝑢^𝑛=𝑢^𝑛ℎ−𝑃_ℎ𝑢^𝑛+𝑃_ℎ𝑢^𝑛−𝑢^𝑛=𝜃^𝑛+𝜌^𝑛,(3.20) where 𝑃_ℎ𝑢^𝑛 is a Rits-Galerkin projection operator defined as follows: 𝑎𝑎(𝑢𝑤)^𝑎𝐷^𝛽𝑥𝑢^𝑛−𝑃_ℎ𝑢^𝑛,𝐷𝑣=0,₀_𝑎𝐷^𝛽𝑥𝑢^𝑛−𝑃_ℎ𝑢^𝑛,𝐷𝑣=0.(3.21)

Lemma 3.5. Let 𝑎(𝑢),𝑏(𝑢) be smooth functions on Ω, 0<𝑚≤𝑎(𝑢),𝑏(𝑢)≤𝑀, and 𝑃_ℎ𝑢^𝑛 is defined as above, then ‖‖_𝑎𝐷‖^𝛼𝑥𝑃𝑢^𝑛_ℎ−𝑃_ℎ𝑢^𝑛‖‖≤𝑐ℎ^{𝑘+1−𝛼}‖𝑢‖_𝑘+1,𝑢^𝑛−𝑢^𝑛‖≤𝑐ℎ^𝑘+1‖𝑢‖_𝑘+1.(3.22)

Proof. Using the definition of 𝑃_ℎ𝑢^𝑛, one gets ‖‖_𝑎𝐷^𝛼𝑥𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖‖²=≤𝑐||‖‖_𝑎𝐷^𝛼𝑥_𝑎𝑃_ℎ𝑢^𝑛−𝑢^𝑛,_𝑎𝐷^𝛼𝑥𝑃_ℎ𝑢^𝑛−𝑢^𝑛||𝐷^𝛼𝑥𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖‖⋅‖‖_𝑎𝐷^𝛼𝑥(𝜒−𝑢^𝑛)‖‖,(3.23) where 𝜒∈𝑆_ℎ. Utilizing the interpolation of 𝐼_ℎ𝑢^𝑛 leads to ‖‖_𝑎𝐷^𝛼𝑥𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖‖≤inf_{𝜒∈𝑆ℎ}𝑐‖𝜒−𝑢‖_𝛼≤𝑐‖‖𝐼_ℎ𝑢^𝑛−𝑢^𝑛‖‖_𝛼≤𝑐ℎ^{𝑘+1−𝛼}‖𝑢‖_𝑘+1.(3.24) Next we estimate ‖𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖. For all𝜙∈𝐿²(Ω), 𝑤 is the solution of the following equation: −_𝑎𝐷𝑤=^𝑥2𝛼𝑤0,=𝜙𝑤,∈𝑤𝜕Ω∈.Ω,(3.25)
So we have ‖𝑤‖_2𝛼≤𝛾₃‖𝜙‖.(3.26) For all𝜒∈𝑆_ℎ, with the help of approximation properties of 𝑆_ℎ and the weak form, we can obtain 𝑃_ℎ𝑢^𝑛−𝑢^𝑛,𝜙==−−𝑃_ℎ𝑥𝑢^𝑛−𝑢^𝑛,_𝑎𝐷^𝑥2𝛼𝑤=−_𝑥𝐷^𝛼𝑏𝑃_ℎ𝑢^𝑛−𝑢^𝑛,_𝑎𝐷^𝛼𝑥𝑤𝐷^𝛼𝑏𝑃_ℎ𝑢^𝑛−𝑢^𝑛,_𝑎𝐷^𝛼𝑥(𝑤−𝜒)≤‖‖𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖‖_𝛼‖𝑤−𝜒‖_𝛼≤‖‖𝑃_ℎ𝑢^𝑛−𝑢^𝑛‖‖_𝛼inf_{𝜒∈𝑆ℎ}‖𝑤−𝜒‖_𝛼‖‖𝑃≤𝑐ℎ^𝑟−𝛼‖𝑢‖_𝑟ℎ^𝛼‖𝑤‖_2𝛼=𝑐ℎ^𝑟‖𝑢‖_𝑟‖𝜙‖,_ℎ𝑢^𝑛−𝑢^𝑛‖‖=sup_{0≠𝜙∈𝐿2(Ω)}𝑃_ℎ𝑢^𝑛−𝑢^𝑛,𝜙‖𝜙‖≤𝑐ℎ^𝑟‖𝑢‖_𝑟.(3.27)

Lemma 3.6 (see [21]). Let 𝑇_ℎ,0<ℎ≤1, denote a quasiuniform family of subdivisions of a polyhedral domain Ω⊂𝑅^𝑑. Let (𝐾^,𝑃,𝑁) be a reference finite element such that 𝑃⊂𝑊^𝑙,𝑝(𝐾^)∩𝑊^𝑚,𝑞(𝐾^) is a finite-dimensional space of functions on 𝐾^,𝑁 is a basis for 𝑃′, where 1≤𝑝≤∞,1≤𝑝≤∞, and 0≤𝑚≤𝑙. For 𝐾∈𝑇_ℎ, let (𝐾,𝑃_𝐾,𝑁_𝐾) be the affine equivalent element, and 𝑉_ℎ=𝑣∶𝑣 is measurable and 𝑣|_𝐾∈𝑃_𝐾,forall𝐾∈𝑇_ℎ. Then there exists a constant 𝐶=𝐶(𝑙,𝑝,𝑞) such that ⎡⎢⎢⎣_{𝑘∈𝑇ℎ}‖𝑣‖^2𝑊_{𝑙,𝑝(𝐾)}⎤⎥⎥⎦^1/𝑝≤𝐶ℎ^{𝑚−𝑙+min(0,𝑑/𝑝−𝑑/𝑞)}⋅⎡⎢⎢⎣_{𝑘∈𝑇ℎ}‖𝑣‖^𝑞𝑊_{𝑚,𝑞(𝐾)}⎤⎥⎥⎦^1/𝑞.(3.28) The following Gronwall’s lemma is useful for the error analysis later on.

Lemma 3.7 (see [2]). Let Δ𝑡,𝐻 and 𝑎_𝑛,𝑏_𝑛,𝑐_𝑛,𝛾_𝑛 (for integer 𝑛≥0) be nonnegative numbers such that 𝑎_𝑁+Δ𝑡^𝑁_𝑛=0𝑏_𝑛≤Δ𝑡^𝑁_𝑛=0𝛾_𝑛𝑎_𝑛+Δ𝑡^𝑁_𝑛=0𝑐_𝑛+𝐻,(3.29) for 𝑁≥0. Suppose that Δ𝑡𝛾_𝑛<1, for all 𝑛, and set 𝜎_𝑛=(1−Δ𝑡𝛾_𝑛)⁻¹. Then 𝑎_𝑁+Δ𝑡^𝑁_𝑛=0𝑏_𝑛≤expΔ𝑡^𝑁_𝑛=0𝜎_𝑛𝛾_𝑛Δ𝑡^𝑁_𝑛=0𝑐_𝑛+𝐻,(3.30) for 𝑁≥0.
The following norms are also used in the analysis: ‖|‖𝑣||‖𝑣|‖_{∞,𝑘0,𝑘}=max_{0≤𝑛≤𝑁}‖𝑣^𝑛‖_𝑘,=^𝑁_𝑛=0𝜏‖𝑣^𝑛‖^2𝑘^1/2.(3.31)

Theorem 3.8. Assume that Problem 1 has a solution 𝑢 satisfying 𝑢_𝑡𝑡,𝑢_𝑡𝑡𝑡∈𝐿²(0,𝑇,𝐿²(Ω)) with 𝑢,𝑢_𝑡∈𝐿²(0,𝑇,𝐻^𝑘+1). If Δ𝑡≤𝑐ℎ, then the finite element approximation is convergent to the solution of Problem 1 on the interval (0,T], as Δ𝑡,ℎ→0. The approximation 𝑢_ℎ also satisfies the following error estimates ‖‖𝑢−𝑢_ℎ‖‖_0,𝛼≤𝐶ℎ^𝑘+1‖‖𝑢_𝑡‖‖_0,𝑘+1+ℎ^{𝑘+1−𝛼}‖𝑢‖_0,𝑘+1+𝜏²‖‖𝑢_𝑡𝑡‖‖_0,0‖‖𝑢−𝑢_ℎ+𝜏ℎ^{𝑘+1−𝛼}‖‖𝑢_𝑡𝑡‖‖_0,𝑘+1+𝜏²‖‖𝑢_𝑡𝑡𝑡‖‖_0,0,(3.32)‖‖_∞,0≤𝐶ℎ^𝑘+1‖‖𝑢_𝑡‖‖_0,𝑘+1+ℎ^{𝑘+1−𝛼}‖𝑢‖_0,𝑘+1+𝜏²‖‖𝑢_𝑡𝑡𝑡‖‖_0,0+𝜏ℎ^{𝑘+1−𝛼}‖‖𝑢_𝑡𝑡‖‖_0,𝑘+1+𝜏²‖‖𝑢_𝑡𝑡‖‖_0,0+ℎ^𝑘+1‖𝑢‖^2∞,𝑘+1.(3.33)

Proof. For 𝑡=𝑡_𝑛−𝜏/2=𝑡_𝑛−1/2,𝑛=0,1,…,𝑁, find 𝑢^𝑛−1/2 such that 𝜕_𝑡𝑢^𝑛−1/2,𝑣+𝑎𝑢^𝑛−1/2_𝑎𝐷_𝑥𝑢^𝑛−1/2,𝐷𝑣−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛−1/2,𝑣=𝑓𝑢^𝑛−1/2,𝑣.(3.34)
Subtracting the above equation from the fully discrete scheme (3.13), and substituting 𝑢^𝑛ℎ−𝑢^𝑛=(𝑢^𝑛ℎ−𝑃_ℎ𝑢^𝑛)+(𝑃_ℎ𝑢^𝑛−𝑢^𝑛)=𝜃^𝑛+𝜌^𝑛 into it, we obtain the following error formulation relating to 𝜃^𝑛 and 𝜌^𝑛: 𝜕_𝑡𝜃=^𝑛,𝑣𝑎+𝑢_𝑛ℎ̃𝑎𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝜃^𝑛,𝐷𝑣−𝑏𝑢^𝑛ℎ̃𝐷𝜃^𝑛,𝑣+𝑎𝑢_𝑎𝐷^𝛽𝑥𝐼_ℎ𝑢^𝑛,𝐷𝑣+𝑏𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝐼_ℎ𝑢^𝑛,𝑣+𝜕_𝑡𝑢^𝑛−1/2,𝑣−𝜕_𝑡𝐼_ℎ𝑢^𝑛,𝑣_𝑛−1/2=−𝑎𝑢_𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2,𝐷𝑣−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛−1/2,𝑣+𝑓𝑢^𝑛ℎ̃,𝑣−𝑓𝑢^𝑛−1/2,𝑣+𝑏𝑢_𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝜌^𝑛,𝐷𝑣+𝑎𝑢^𝑛−1/2_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2−𝑎𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝑢^𝑛,𝐷𝑣+𝑓𝑢_𝑛ℎ̃𝐷𝜌^𝑛,𝑣+𝑏𝑢^𝑛ℎ̃𝐷𝑢^𝑛−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛−1/2,𝐷𝑣=𝑅₁−𝑓𝑢^𝑛−1/2,𝑣+𝜕_𝑡𝑢^𝑛−1/2−𝜕_𝑡𝐼_ℎ𝑢^𝑛,𝑣(𝑣)+𝑅₂(𝑣)+𝑅₃(𝑣)+𝑅₄(𝑣)+𝑅₅(𝑣)+𝑅₆(𝑣).(3.35) Setting 𝑣=𝜃^𝑛, we obtain 𝜕_𝑡𝜃=^𝑛,𝑅𝜃^𝑛₁+𝑎(𝑢^𝑛̃)_𝑎𝐷^𝛽𝑥𝜃^𝑛,𝐷𝜃^𝑛−𝑏(𝑢^𝑛̃)𝐷𝜃^𝑛,𝜃^𝑛𝜃^𝑛+𝑅²𝜃^𝑛+𝑅³𝜃^𝑛+𝑅⁴𝜃^𝑛+𝑅⁵𝜃^𝑛+𝑅⁶𝜃^𝑛.(3.36) Note that 𝜕_𝑡𝜃^𝑛,𝜃^𝑛=𝜃^𝑛−𝜃^𝑛−1𝜏,𝜃^𝑛+𝜃^𝑛−12=21𝜏‖𝜃^𝑛‖²−‖‖𝜃^𝑛−1‖‖².(3.37) According to (3.2) and Lemma 3.2, we have 𝑎(𝑢^𝑛̃)_𝑎𝐷^𝛽𝑥𝜃^𝑛,𝐷𝜃^𝑛≥𝑚|||𝜃^𝑛|||^2𝛼≥𝑐||𝜃^𝑛||^2𝛼+||𝜃^𝑛−1||^2𝛼.(3.38) From Lemma 3.3, the following inequality can be derived: 𝑏(𝑢^𝑛̃)𝜃^𝑛,𝐷𝜃^𝑛≤𝑐⁰𝜀^−𝑐1₂‖𝐷𝑏(𝑢^𝑛̃)‖^𝑐2‖‖𝜃^𝑛‖‖²+𝜀₃‖‖𝜃^𝑛‖‖^2𝛼=𝑐₀𝜀^−𝑐1₂‖𝐷𝑏(𝑢^𝑛̃)‖^𝑐2‖‖‖𝜃^𝑛+𝜃^𝑛−12‖‖‖²+𝜀₃‖‖‖𝜃^𝑛+𝜃^𝑛−12‖‖‖^2𝛼≤𝑐₃𝜀^−𝑐1₂‖𝐷𝑏(𝑢^𝑛̃)‖^𝑐2‖𝜃^𝑛‖²+‖‖𝜃^𝑛−1‖‖²+𝑐⁴𝜀₃‖𝜃^𝑛‖^2𝛼+‖‖𝜃^𝑛−1‖‖^2𝛼.(3.39)
Substituting (3.37)–(3.39) into (3.36) then multiplying (3.36) by 2𝜏, summing from 𝑛=1 to 𝑁, we have ‖𝜃^𝑛‖²−‖‖𝜃²‖‖²+𝜏^𝑁_𝑛=12𝑚𝑐−2𝑐₄≤2𝜏_𝑁𝜀₃‖𝜃^𝑛‖^2𝛼+‖‖𝜃^𝑛−1‖‖^2𝛼_𝑛=1+2𝜏_𝑁𝑐₃𝜀^−𝑐1₂‖‖𝐷𝑏𝑢^𝑛ℎ̃‖‖^𝑐2‖𝜃^𝑛‖²+‖‖𝜃^𝑛−1‖‖²_𝑛=3𝑅¹𝜃^𝑛+𝑅²𝜃^𝑛+𝑅³𝜃^𝑛+𝑅⁴𝜃^𝑛+𝑅⁵𝜃^𝑛+𝑅⁶𝜃^𝑛.(3.40) We now estimate 𝑅₁ to 𝑅₆ in the right hand of (3.40), 𝑅₁𝜃^𝑛≤=𝑀^𝑎𝐷^𝑥1−𝛼_𝑎𝑎𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝜌^𝑛,^𝑎𝐷^𝛼𝑥𝜃^𝑛≤𝜀₄𝐷^𝛼𝑥𝜌^𝑛,_𝑎𝐷^𝛼𝑥𝜃^𝑛≤𝑀‖‖^𝑎𝐷^𝛼𝑥𝜌^𝑛‖‖‖‖_𝑎𝐷^𝛼𝑥𝜃^𝑛‖‖‖‖𝜃^𝑛‖‖^2𝛼+𝑐²⁵4𝜀₄‖‖𝜌^𝑛‖‖^2𝛼=𝜀⁴2‖‖𝜌^𝑛+𝜌^𝑛−1‖‖^2𝛼+𝑐²⁵16𝜀₄‖‖𝜃^𝑛+𝜃^𝑛−1‖‖^2𝛼≤𝜀₄𝑐₆‖𝜃^𝑛‖^2𝛼+‖‖𝜃^𝑛−1‖‖^2𝛼+𝑐⁷𝜀₄‖𝜌^𝑛‖^2𝛼+‖‖𝜌^𝑛−1‖‖^2𝛼.(3.41)
Secondly, we deduce the estimation of 𝑅₂, 𝑅₂𝜃^𝑛==−𝑎𝑎𝑢𝑢^𝑛ℎ̃_{𝑎𝑛−1/2}𝐷^𝛽𝑥𝑢^𝑛,𝐷𝜃^𝑛+𝑎𝑢^𝑛−1/2_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2,𝐷𝜃^𝑛=𝑅₂₁−𝑎𝑢^𝑛ℎ̃_𝑎𝐷^𝛽𝑥𝑢^𝑛,𝐷𝜃^𝑛+𝑎𝑢^𝑛−1/2_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2−_𝑎𝐷^𝛽𝑥𝑢^𝑛,𝐷𝜃^𝑛+𝑅₂₂,(3.42) where 𝑅₂₁≤=𝑐₈𝑎𝑢^𝑛−1/2−𝑎(𝑢^𝑛̃)_𝑎𝐷^𝛽𝑥𝑢^𝑛̃,𝐷𝜃^𝑛4𝜀₅‖‖𝑎𝑢^𝑛−1/2−𝑎(𝑢^𝑛̃)_𝑎𝐷^𝛽𝑥𝑢^𝑛̃‖‖^21−𝛼+𝜀₅‖‖𝐷𝜃^𝑛‖‖^2𝛼−1≤𝑐₉‖‖𝑎𝑢^𝑛−1/2−𝑎(𝑢^𝑛̃)‖‖‖‖_𝑎𝐷^𝛽𝑥𝑢^𝑛̃‖‖^21−𝛼+𝜀₅‖‖𝜃^𝑛‖‖^2𝛼𝑅₂₂≤𝑐₉𝐿‖‖𝑢^𝑛−1/2−𝑢^𝑛̃‖‖+𝜀₅‖‖𝜃^𝑛‖‖^2𝛼,=≤𝑐𝑎₁₀𝑢^𝑛−1/2_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2−_𝑎𝐷^𝛽𝑥𝑢^𝑛,𝐷𝜃^𝑛4𝜀₆‖‖𝑎𝑢^𝑛−1/2_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2−_𝑎𝐷^𝛽𝑥𝑢^𝑛‖‖^21−𝛼+𝜀₆‖‖𝐷𝜃^𝑛‖‖^2−1+𝛼≤𝑐₁₀‖‖𝑎𝑢^𝑛−1/2‖‖²‖‖_𝑎𝐷^𝛽𝑥𝑢^𝑛−1/2−_𝑎𝐷^𝛽𝑥𝑢^𝑛‖‖^21−𝛼+𝜀₆‖‖𝜃^𝑛‖‖^2𝛼≤𝑐₁₀𝑀²‖‖𝑢^𝑛−1/2−𝑢^𝑛‖‖^2𝛼+𝜀₆𝑐‖𝜃^𝑛‖^2𝛼+‖‖𝜃^𝑛−1‖‖^2𝛼.(3.43) The estimations of ‖𝑢^𝑛̃−𝑢^𝑛−1/2‖ and ‖𝑢^𝑛−𝑢^𝑛−1/2‖_𝛼 can be derived as follows: ‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖=‖‖‖‖32𝑢^𝑛−1/2−𝜏2𝑢^{𝑡𝑛−1/2}+𝑢^{𝑛−1/2𝑡𝑡}2!𝜏2²−1+𝑂𝜏³2𝑢^𝑛−1/2−3𝜏2𝑢^{𝑡𝑛−1/2}+𝑢^{𝑛−1/2𝑡𝑡}2!32𝜏²+𝑂𝜏³−𝑢^𝑛−1/2≤𝑐₁₁‖‖‖‖𝜏²‖‖𝑢_𝑡𝑡𝑡_𝑛−1/2‖‖≤𝑐₁₁𝜏²^𝑡𝑛_𝑡𝑛−1‖‖‖‖𝑢_𝑡𝑡(⋅,𝑠)‖‖𝑑𝑠,𝑢^𝑛−𝑢^𝑛−1/2‖‖_𝛼=‖‖‖𝜏⁻¹^𝑡𝑛_{𝑡𝑛−1/2}𝑠−𝑡_𝑛²𝑢_𝑡𝑡(𝑠)𝑑𝑠+^{𝑡𝑛−1/2}_𝑡𝑛−1𝑠−𝑡_𝑛−1²𝑢_𝑡𝑡(𝑠)𝑑𝑠‖‖‖_𝛼≤𝑐₁₂𝜏‖‖‖^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡(𝑠)𝑑𝑠‖‖‖_𝛼≤≤𝑐𝑐₁₂₁₂𝜏^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡𝑡(𝑠)‖‖_𝛼𝑑𝑠𝜏ℎ^{𝑘+1−𝛼}^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡𝑡(𝑠)‖‖_𝑘+1𝑑𝑠.(3.44) Thirdly, it is turn to consider 𝑅₃, 𝑅₃𝜃^𝑛=𝑏𝑢^𝑛ℎ̃𝐷𝜌^𝑛,𝜃^𝑛≤‖‖𝑏𝑢^𝑛ℎ̃𝐷𝜌^𝑛‖‖_−𝛼‖‖𝜃^𝑛‖‖_𝛼≤𝑐¹³4𝜀₇‖‖𝑏𝑢^𝑛ℎ̃‖‖²‖‖𝜌^𝑛‖‖^21−𝛼+𝜀₇‖‖𝜃^𝑛‖‖^2𝛼≤𝑐¹⁴4𝜀₇‖𝜌^𝑛‖^21−𝛼+‖‖𝜌^𝑛−1‖‖^21−𝛼+𝜀⁷𝑐₁₅‖𝜃^𝑛‖^2𝛼+‖‖𝜃^𝑛−1‖‖^2𝛼.(3.45) Next, 𝑅₄𝜃^𝑛==𝑏𝑏𝑢𝑢^𝑛ℎ̃_𝑛ℎ̃𝐷𝑢^𝑛,𝜃^𝑛−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛−1/2,𝜃^𝑛=𝑅₄₁𝐷𝑢^𝑛−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛,𝜃^𝑛+𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛−𝐷𝑢^𝑛−1/2,𝜃^𝑛+𝑅₄₂,(3.46) where 𝑅₄₁≤𝑐¹⁶4𝜀₈‖‖𝑏𝑢^𝑛ℎ̃−𝑏𝑢^𝑛−1/2𝐷𝑢^𝑛‖‖^21−𝛼+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼≤𝑐¹⁶𝐿4𝜀₈‖‖𝑢^𝑛ℎ̃−𝑢^𝑛−1/2‖‖²‖‖𝐷𝑢^𝑛‖‖^21−𝛼+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼=𝑐¹⁶𝐿4𝜀₈‖‖𝑢^𝑛ℎ̃−𝑢^𝑛̃+𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²‖‖𝐷𝑢^𝑛‖‖^21−𝛼+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼≤𝑐₁₇‖‖𝑢^𝑛ℎ̃−𝑢^𝑛̃‖‖²+𝑐₁₇‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼≤𝑐₁₇‖‖𝜃̃^𝑛+𝜌̃^𝑛‖‖²+𝑐₁₇‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼≤𝑐₁₈‖‖𝜃̃^𝑛‖‖²+‖𝜌̃^𝑛‖²+𝑐¹⁷‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²+𝜀₈‖‖𝜃^𝑛‖‖^2𝛼.(3.47) Rewriting 𝑅₄₂ by the aid of (3.20), we have 𝑅₄₂≤𝑐¹⁹4𝜀₉‖‖𝑢^𝑛−𝑢^𝑛−1/2‖‖²+𝜀₉‖‖𝜃^𝑛‖‖^2𝛼.(3.48) The estimation of 𝑅₅ is deduced as follows: 𝑅₅𝜃^𝑛≤≤𝐿‖‖𝑓‖‖𝑢_𝑛ℎ̃𝑢^𝑛ℎ̃−𝑓𝑢^𝑛−1/2‖‖‖‖𝜃^𝑛‖‖≤𝐿𝑐₂₀−𝑢^𝑛−1/2‖‖‖‖𝜃^𝑛‖‖4𝜀₁₀‖‖𝑢^𝑛ℎ̃−𝑢^𝑛−1/2‖‖²+𝜀₁₀‖‖𝜃^𝑛‖‖²≤𝐿𝑐₂₁‖‖𝜃̃^𝑛+𝜌̃^𝑛‖‖²+‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²+𝜀¹⁰‖‖𝜃^𝑛‖‖²≤𝑐₂₂‖‖𝜃̃^𝑛‖‖²+‖𝜌̃^𝑛‖²+𝐿𝑐²¹‖‖𝑢^𝑛̃−𝑢^𝑛−1/2‖‖²+𝜀₁₀‖‖𝜃^𝑛‖‖².(3.49) Last, we estimate 𝑅₆, 𝑅₆𝜃^𝑛==𝜕𝜕^𝑡_𝑡𝑢^𝑛−1/2,𝜃^𝑛−𝜕_𝑡𝑃_ℎ𝑢^𝑛,𝜃^𝑛=𝜕_𝑡𝑢^𝑛−1/2−𝜕𝑢^𝑛,𝜃^𝑛+𝜕_𝑡𝑢^𝑛−𝜕_𝑡𝑃_ℎ𝑢^𝑛,𝜃^𝑛≤‖‖𝜕_𝑡𝑢^𝑛−1/2−𝜕_𝑡𝑢^𝑛,𝜃^𝑛+𝜕_𝑡𝜌^𝑛,𝜃^𝑛𝑢^𝑛−1/2−𝜕_𝑡𝑢^𝑛‖‖‖𝜃^𝑛‖+‖‖𝜕_𝑡𝜌^𝑛‖‖‖𝜃^𝑛‖,(3.50) where ‖‖𝜕_𝑡𝑢^𝑛−1/2−𝜕𝑢^𝑛‖‖=(2𝜏)⁻¹𝑐₂₃‖‖‖^𝑡𝑛_{𝑡𝑛−1/2}≤𝑐₂₃𝑠−𝑡_𝑛²𝑢_𝑡𝑡𝑡(𝑠)𝑑𝑠+^{𝑡𝑛−1/2}_𝑡𝑛−1𝑠−𝑡_𝑛−1²𝑢_𝑡𝑡𝑡(𝑠)𝑑𝑠‖‖‖≤𝑐₂₃𝜏‖‖‖^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡𝑡(𝑠)𝑑𝑠‖‖‖‖‖𝜏^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡𝑡𝑡(𝑠)‖‖𝑑𝑠,𝜕_𝑡𝜌^𝑛‖‖=‖‖‖𝜌^𝑛−𝜌^𝑛−1𝜏‖‖‖≤𝜏⁻¹‖‖‖^𝑡𝑛_𝑡𝑛−1≤𝜏₋₁𝜌^𝑛𝑡(𝑠)𝑑𝑠‖‖‖^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡(𝑠)‖‖𝑑𝑠≤𝜏⁻¹^𝑡𝑛_𝑡𝑛−1=1𝑑𝑠^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡(𝑠)‖‖𝑑𝑠_{𝑡𝑛𝑡𝑛−1}‖‖𝑢_𝑡(𝑠)‖‖≤ℎ^𝑘+1^𝑡𝑛_𝑡𝑛−1‖‖𝑢_𝑡(𝑠)‖‖_𝑘+1𝑑𝑠.(3.51)
The ‖𝜃²‖ should be estimated with (3.14). Let 𝑛=1 then subtracting (3.34) from the two equations of (3.14), respectively, one gets 𝜕_𝑡𝜃=^1,0−,𝑣𝑎+𝑢_0ℎ𝑎𝑢^0ℎ_𝑎𝐷^𝛽𝑥𝜃^1,0,𝐷𝑣−𝑏𝑢^0ℎ𝐷𝜃^1,0,𝑣+𝑏𝑢_𝑎𝐷^𝛽𝑥𝜌^1,0,𝐷𝑣−𝑎𝑢^1/2_𝑎𝐷^𝛽𝑥𝑢^1/2−𝑎𝑢^0ℎ_𝑎𝐷^𝛽𝑥𝑢^1,0,𝐷𝑣_0ℎ+𝑓𝑢_0ℎ𝐷𝜌^1,0,𝑣+𝑏𝑢^0ℎ𝐷𝑢^1,0−𝑏𝑢^1/2𝐷𝑢^1/2,𝐷𝑣=𝑅₁−𝑓𝑢^1/2,𝑣+𝜕_𝑡𝑢^1/2−𝜕_𝑡𝐼_ℎ𝑢^1,0,𝑣(𝑣)+𝑅₂(𝑣)+𝑅₃(𝑣)+𝑅₄(𝑣)+𝑅₅(𝑣)+𝑅₆(𝑣),𝜕_𝑡𝜃=¹,−𝑣+𝑎𝑎𝑢_0ℎ𝑢^0ℎ_𝑎𝐷^𝛽𝑥𝜃¹,𝐷𝑣−𝑏𝑢^0ℎ𝐷𝜃¹,𝑣+𝑢^ℎ1,02^𝑎𝐷^𝛽𝑥𝜌¹,𝐷𝑣−𝑎𝑢^1/2_𝑎𝐷^𝛽𝑥𝑢^1/2−𝑎𝑢^0ℎ+𝑢^ℎ1,02^𝑎+𝑏𝑢𝐷^𝛽𝑥𝑢¹,𝐷𝑣_0ℎ+𝑢^ℎ1,02𝐷𝜌¹,𝑣+𝑏𝑢^0ℎ+𝑢^ℎ1,02𝐷+𝑓𝑢_0ℎ𝑢¹−𝑏𝑢^1/2𝐷𝑢^1/2,𝐷𝑣+𝑢^ℎ1,02−𝑓𝑢^1/2,𝑣+𝜕_𝑡=𝑅₁𝑢^1/2−𝜕_𝑡𝐼_ℎ𝑢¹,𝑣(𝑣)+𝑅₂(𝑣)+𝑅₃(𝑣)+𝑅₄(𝑣)+𝑅₅(𝑣)+𝑅₆(𝑣).(3.52)
Setting 𝑣=𝜃^1,0, and using the similar estimation (see (3.40)), one has ‖‖𝜃^1,0‖‖²≤𝑐𝜏²^𝑡1_𝑡0‖‖𝑢_𝑡𝑡(𝑠)‖‖^2𝛼𝑑𝑠+ℎ^{2(𝑘+1−𝛼)}‖𝑢‖^2𝑘+1+𝜏²^𝑡1_𝑡0+𝑐ℎ_2(𝑘+1)‖‖𝑢_𝑡𝑡𝑡(𝑠)‖‖²𝑑𝑠^𝑡1_𝑡0‖‖𝑢_𝑡(𝑠)‖‖^2𝑘+1𝑑𝑠.(3.53)
Letting 𝑣=𝜃¹, applying the above result of 𝜃^1,0, and using the similar estimation (see (3.53)), we get ‖‖𝜃¹‖‖²≤𝑐𝜏²^𝑡1_𝑡0‖‖𝑢_𝑡𝑡(𝑠)‖‖^2𝛼𝑑𝑠+ℎ^{2(𝑘+1−𝛼)}‖𝑢‖^2𝑘+1+𝜏²^𝑡1_𝑡0+𝑐ℎ_2(𝑘+1)‖‖𝑢_𝑡𝑡𝑡(𝑠)‖‖²𝑑𝑠^𝑡1_𝑡0‖‖𝑢_𝑡(𝑠)‖‖^2𝑘+1𝑑𝑠.(3.54)
Using 𝑇=𝑁𝜏 and Gronwall’s lemma, we get ‖‖||𝜃||‖‖^20,𝛼=^𝑁_𝑛=0𝜏‖𝜃‖^2𝛼.(3.55)
Hence, using the interpolation property and ‖‖||𝑢−𝑢_ℎ||‖‖_0,𝛼≤‖‖||𝜃||‖‖_0,𝛼+‖‖||𝜌||‖‖_0,𝛼,(3.56) the estimate (3.32) holds.
Also using the interpolation property, Gronwall’s lemma, and the approximation properties, we get ‖‖||𝑢−𝑢_ℎ||‖‖_∞,0≤‖‖||𝜃||‖‖_∞,0+‖‖||𝜌||‖‖_∞,0≤max_{0≤𝑛≤𝑁}‖‖||𝜃^𝑛||‖‖²+ℎ^2𝑘+1‖|𝑢|‖^2∞,𝑘+1,(3.57) which is just the estimate (3.33).

4. Numerical Examples

In this section, we present the numerical results which confirm the theoretical analysis in Section 3.

Let 𝐾 denote a uniform partition on [0,𝑎], and 𝑆_ℎ the space of continuous piecewise linear functions on 𝐾, that is, 𝑘=1.. In order to implement the Galerkin finite element approximation, we adapt finite element discrete along the space axis, and finite difference scheme along the time axis. We associate shape function of space 𝑋_ℎ with the standard basis of hat functions on the uniform grid of size ℎ=1/𝑛. We have the predicted rates of convergence if the condition Δ𝑡=𝑐ℎ of ‖‖‖‖𝑢𝑢−−𝑢𝑢_ℎℎ‖‖_0,𝛼∼𝑂ℎ^2−𝛼,‖‖_∞,0∼𝑂ℎ^2−𝛼,(4.1) provided that the initial value 𝜑(𝑥) is smooth enough.

Example 4.1. The following equation 𝜕𝑢𝜕𝑡=𝐷𝑢²⁰𝐷^𝑥0.5𝑢(𝑥,𝑡)−2𝑥(𝑥−1)2𝑥^1.5Γ(2.5)−Γ𝑥(1.^0.55)𝑒^−2𝑡−𝑢₂𝐷𝑢−𝑢(𝑥,𝑡)𝑒^−𝑡Γ2(𝑥1.5𝑢^0.5𝑢)((0𝑥−,,𝑡0Γ)𝑥)=(=0.𝑢5^−0.5𝑥)(1(𝑥,𝑡,−)=10)0≤,,𝑥00≤≤≤1𝑥,𝑡≤≤0≤11,,𝑡≤1,(4.2) has a unique solution 𝑢(𝑥,𝑡)=𝑒^−𝑡𝑥(𝑥−1).
If we select Δ𝑡=𝑐ℎ and note that the initial value 𝑢⁰ is smooth enough, then we have ‖‖‖‖𝑢𝑢−−𝑢𝑢_ℎℎ‖‖_0,0.75∼𝑂ℎ^1.25,‖‖_∞,0∼𝑂ℎ^1.25.(4.3)

Table 1 includes numerical calculations over a regular partition of [0,1]. We can observe the experimental rates of convergence agree with the theoretical rates for the numerical solution.

Table 1 

Numerical error result for Example 4.1.

Example 4.2. The function 𝑢(𝑥,𝑡)=cos(𝑡)𝑥²(2−𝑥)² solves the equation in the following form: 𝜕𝑢𝜕𝑡=⁰𝐷^𝑥1.7𝑢(𝑥,𝑡)+𝑏(𝑢)𝑢𝐷(𝑢𝑥−,0𝑢)[=4𝑥(1₂−𝑥)+tan𝑡]+𝑓(𝑥,𝑡),𝑥∈(0,2),𝑡∈[0,1),𝑢(0,𝑡)=0,𝑢(2(2−,𝑡𝑥))=²,0,0≤0≤𝑥≤𝑡≤2,1,(4.4) where 𝑏(𝑢)=√𝑢𝑓(𝑥,𝑡)=cos𝑡√cos𝑡,cos(0.85𝜋)24𝑥^2.3+(2−𝑥)^2.3Γ(3.3)−24𝑥^1.3+(2−𝑥)^1.3Γ(2.3)−8𝑥^0.3+(2−𝑥)^0.3Γ(1.3).(4.5)
If we select Δ𝑡=𝑐ℎ, then ‖‖‖‖𝑢𝑢−−𝑢𝑢_ℎℎ‖‖_0,0.85∼𝑂ℎ^1.15,‖‖_∞,0∼𝑂ℎ^1.15.(4.6)

Table 2 shows the error results at different size of space grid. We can observe that the experimental rates of convergence still support the theoretical rates.

Table 2 

Numerical error result for Example 4.2.

Example 4.3. Consider the following space-fractional differential equation with the nonhomogeneous boundary conditions, 𝜕𝑢𝜕𝑡=⁰𝐷^𝑥1.7𝑢(𝑥,𝑡)−3𝑥^𝑥0𝑢𝑑𝑥−2𝑥^0.3𝑒^−𝑡𝑢(𝑥,0)=𝑥₂Γ(1.3),0≤𝑥≤1,0≤𝑡≤1,𝑢(0,𝑡)=0,𝑢(1,𝑡,)=0𝑒≤_−𝑡𝑥≤1,,0≤𝑡≤1,(4.7) whose exact solution is 𝑢(𝑥,𝑡)=𝑒^−𝑡𝑥².
We still choose Δ𝑡=𝑐ℎ, then get the convergence rates ‖‖‖‖𝑢𝑢−−𝑢𝑢_ℎℎ‖‖_0,0.85∼𝑂ℎ^1.15,‖‖_∞,0∼𝑂ℎ^1.15.(4.8)

The numerical results are presented in Table 3 which are in line with the theoretical analysis.

Table 3 

Numerical error result for Example 4.3.

5. Conclusion

In this paper, we propose a fully discrete Galerkin finite element method to solve a type of fractional advection-diffusion equation numerically. In the temporal direction we use the modified Crank-Nicolson method, and in the spatial direction we use the finite element method. The error analysis is derived on the basis of fractional derivative space. The numerical results agree with the theoretical error estimates, demonstrating that our algorithm is feasible.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China under grant no. 10872119, the Key Disciplines of Shanghai Municipality under grant no. S30104, the Key Program of Shanghai Municipal Education Commission under grant no. 12ZZ084, and the Natural Science Foundation of Anhui province KJ2010B442.