A Finite Element Approximation for a Class of Caputo Time-Fractional Diffusion Equations

TL;DR

This paper introduces a finite element and finite difference scheme for solving Caputo time-fractional diffusion equations, providing stability, error analysis, and numerical validation of the method's accuracy and efficiency.

Contribution

It presents a novel fully discrete numerical scheme combining finite difference in time and finite element in space for Caputo fractional diffusion equations, with proven stability and error estimates.

Findings

The method is stable under certain conditions.

Error estimates demonstrate the scheme's accuracy.

Numerical examples confirm efficiency and reliability.

Abstract

We develop a fully discrete scheme for time-fractional diffusion equations by using a finite difference method in time and a finite element method in space. The fractional derivatives are used in Caputo sense. Stability and error estimates are derived. The accuracy and efficiency of the presented method is shown by conducting two numerical examples.

Tables

Table 1. Table 1. Results for the problem of Example  4.1 with α = 0.1 , 0.5 , 0.9 𝛼 0.1 0.5 0.9 \alpha=0.1,0.5,0.9 .

$x_i$	exact sol	$\alpha$	approximate sol	error
$x_{100}$		0.1	$0.571093$	${2.10}^{-6}$
	$0.571095$	0.5	$0.571100$	$-{5.10}^{-6}$
		0.9	$0.571175$	$-{8.10}^{-5}$
$x_{300}$		0.1	$0.934012$	${3.10}^{-6}$
	0.934015	0.5	$0.934023$	$-{8.10}^{-6}$
		0.9	$0.934145$	$-1,{3.10}^{-4}$
$x_{500}$		0.1	$0.006158$	$0$
	0.006158	0.5	$0.006158$	$0$
		0.9	$0.006159$	$-{1.10}^{-6}$
$x_{700}$		0.1	$-0.930207$	$-{3.10}^{-6}$
	-0.930209	0.5	$-0.930217$	${8.10}^{-6}$
		0.9	$-0.930339$	$1,{3.10}^{-4}$
$x_{900}$		0.1	$-0.581057$	$-{2.10}^{-6}$
	-0.581059	0.5	$-0.581064$	${5.10}^{-6}$
		0.9	$-0.581140$	$8,{1.10}^{-6}$

Equations

$$\begin{array}{rl}& {}_0^C{D}_t^{\alpha }u(x,t)-\Delta u(x,t)=f(x,t),x\in \Omega ,t\in [0,T],\\ & u(x,0)=u_0(x),x\in \Omega ,\\ & u(x,t)=0,x\in \partial \Omega ,t\in [0,T],\end{array}$$

$${}_0^C{D}_t^{\alpha }u(x,t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial u(x,s)}{\partial s}\frac{ds}{(t-s{)}^{\alpha }},0<\alpha <1,$$

$${}_0^C{D}_t^{\alpha }u(t)+Au(t)=f(t),0<t<T,0<\alpha <1,$$

$$u(0)=u_0.$$

$$\begin{array}{rl}{}_0^C{D}_t^{\alpha }u(t_{k+1})& =\frac{1}{\Gamma (2-\alpha )\Delta t^{\alpha }}\sum _{j=0}^k\left(u(t_{k+1-j})-u(t_{k-j})\right)((j+1{)}^{1-\alpha }-j^{1-\alpha })+{\tilde{R}}_{k+1},\end{array}$$

$${\tilde{R}}_{k+1}\le c_u\Delta t^{2-\alpha }$$

$$\begin{array}[]{lll}&b_{j}>0,\quad j=0,1,\ldots,k,\\ &1=b_{0}>b_{1}>\cdots>b_{k},\quad b_{k}\longrightarrow 0\quad as\leavevmode\nobreak\ k\longrightarrow\infty,\\ &\sum_{j=0}^{k}(b_{j}-b_{j+1})+b_{k+1}=(1-b_{1})+\sum_{j=1}^{k-1}(b_{j}-b_{j+1})+b_{k}=1.\end{array}$$

$$L_t^{\alpha }u(t_{k+1})=\frac{1}{{\alpha }_0}\left(u(t_{k+1})-(1-b_1)u(t_k)-\sum _{j=1}^{k-1}(b_j-b_{j+1})u(t_{k-j})-b_{k}u(t_0)\right)$$

$${\alpha }_0:=\Gamma (2-\alpha )\Delta t^{\alpha }.$$

$$\begin{array}{rl}{}_0^C{D}_t^{\alpha }u(t_{k+1})& =\frac{1}{{\alpha }_0}\left(u(t_{k+1})-(1-b_1)u(t_k)-\sum _{j=1}^{k-1}(b_j-b_{j+1})u(t_{k-j})-b_{k}u(t_0)\right)+{\tilde{R}}_{k+1}\\ & =L_t^{\alpha }u(t_{k+1})+{\tilde{R}}_{k+1}.\end{array}$$

$$L_t^{\alpha }u(t_{k+1})+Au(t_{k+1})=f(t_{k+1})-{\tilde{R}}_{k+1},k=0,1,\dots ,K-1.$$

$$L_t^{\alpha }{u}^{k+1}+A{u}^{k+1}=f^{k+1},k=0,1,\dots ,K-1.$$

$$u^{k+1}(x)=0,k\ge 0,x\in \partial \Omega ,$$

$$u^0(x)=u_0(x),x\in \Omega .$$

$$u^{k+1}+{\alpha }_{0}A{u}^{k+1}={\alpha }_0{f}^{k+1}+(1-b_1)u^k+\sum _{j=1}^{k-1}(b_j-b_{j+1})u^{k-j}+b_k{u}^0,k\ge 1,$$

$$u^1+{\alpha }_{0}A{u}^1={\alpha }_0{f}^1+u^0.$$

$$r^{k+1}:={\alpha }_0\left({}_0^C{D}_t^{\alpha }u(t_{k+1})-L_t^{\alpha }u(t_{k+1})\right),$$

$$|r^{k+1}|={\alpha }_0|r_{\Delta t}^{k+1}|\le c_u\Delta t^2.$$

$$H^1(\Omega ):=\left\{v\in L^2(\Omega ),\nabla v\in L^2(\Omega )\right\},$$

$$H_0^1(\Omega ):=\left\{v\in H^1(\Omega ),v{|}_{\partial \Omega }=0\right\},$$

$$H^m(\Omega ):=\left\{v\in L^2(\Omega ),\frac{d^{k}v}{d{x}^k}\in L^2(\Omega )\text{leavevmode}\text{ for all positive integer }k\le m\right\},$$

$$(u,v)={\int }_{\Omega }uvdx,(u,v{)}_1=(u,v)+(\nabla u,\nabla v),$$

$$\Vert v\Vert =(v,v{)}^{\frac{1}{2}},\Vert v{\Vert }_1=(v,v{)}_1^{\frac{1}{2}},\text{and}\Vert v{\Vert }_2=\Vert v{\Vert }_{H^2}={\left(\sum _{k\le 2}{\Vert \frac{d^{k}v}{d{x}^k}\Vert }^2\right)}^{1/2}.$$

$$(u^{k+1},v)+{\alpha }_0(A{u}^{k+1},v)={\alpha }_0(f^{k+1},v)+(1-b_1)(u^k,v)+\sum _{j=1}^{k-1}(b_j-b_{j+1})(u^{k-j},v)+b_k(u^0,v)$$

$$\Vert u^k\Vert \le \Vert u^0\Vert +{\alpha }_{0}k\Vert f{\Vert }_{L^{\infty }},k=1,2,\dots ,K.$$

$$(I+{\alpha }_{0}A)u^1={\alpha }_0{f}^1+u^0.$$

$$u^1=(I+{\alpha }_{0}A{)}^{-1}({\alpha }_0{f}^1+u^0).$$

$$\Vert (I+{\alpha }_{0}A{)}^{-1}\Vert =\underset{{\lambda }_j>0}{\sup }|(1+{\alpha }_0{\lambda }_j{)}^{-1}|<1.$$

$$\Vert u^1\Vert \le \Vert u^0+{\alpha }_0{f}^1\Vert .$$

$$\Vert u^1\Vert \le \Vert u^0\Vert +{\alpha }_0\Vert f{\Vert }_{L^{\infty }},$$

$$\Vert u^j\Vert \le \Vert u^0\Vert +j{\alpha }_0\Vert f{\Vert }_{L^{\infty }},\forall j=1,2,\dots ,k.$$

Full text

∗Corresponding author: Delfim F. M. Torres (delfim@ua.pt)

Moulay Rchid Sidi Ammi and Ismail Jamiai

AMNEA Group, Department of Mathematics,

Faculty of Sciences and Technics, Moulay Ismail University,

B.P. 509, Errachidia, Morocco

Delfim F. M. Torres∗

Center for Research and Development in Mathematics and Applications (CIDMA),

Department of Mathematics, University of Aveiro,

3810-193 Aveiro, Portugal

1. Introduction

Fractional calculus is the field of mathematical analysis that deals with the investigation and application of integrals and derivatives of arbitrary order. The fractional calculus may be considered an old topic, starting from some speculations of Leibniz and Euler, respectively in the 17th and 18th centuries, and yet a recent subject under strong development [2, 3, 22].

In recent years, time-fractional partial differential equations (TFPDEs) have aroused a considerable interest among mathematicians and also have been applied broadly in various applications of numerical analysis in different research areas, including fractal phenomena, diffusion processes, complex networks, stochastic interfaces, synoptic climatology, option pricing mechanisms, medical image processing, electromagnetic, electro-chemistry and material sciences, and chaotic dynamics of nonlinear systems [6, 14, 23]. In view of the importance of TFPDEs, many researchers investigate them in both analytical and numerical frameworks. Several works and methods have been developed, such as finite difference methods [5, 18, 19, 20, 26], finite element methods [7, 9], spectral methods [15], Adomian decomposition methods [21], and variational iteration methods [8]. Regarding analytical solutions to TFDEs, one can use Green and Fox functions and their properties, similarity methods, and Fourier–Laplace transforms or Wright functions [11, 12, 13, 17].

Here we study a numerical approach to the following initial-boundary value time-fractional Caputo diffusion problem:

[TABLE]

where $\alpha$ is the order of the time-fractional derivative, $0<\alpha<1$, and $\Omega$ is a bounded open domain in $\mathbb{R}^{d}$, $1\leq d\leq 3$. The operator ${}^{C}_{0}D^{\alpha}_{t}$ is the Caputo fractional derivative of order $\alpha$ of function $u(x,t)$, defined by

[TABLE]

where $\Gamma$ denotes the Gamma function. In [16], some analytical solutions of the time-fractional diffusion equation (1) with a vanishing forcing term (i.e., $f(x,t)\equiv 0$) are obtained, by applying the finite sine and Laplace transforms based on the fundamental Mittag–Leffler function. It is a hard task to search and to compute the exact solution, especially for large time, due to slow convergence of the series of the Mittag–Leffler function. Therefore, developing efficient numerical methods is a significant question, and considerable efforts have been devoted to develop numerical algorithms for this class of problems. In general, finite difference methods and finite element methods are the most accepted approaches for solving FPDEs. For instance, in [16] a practical finite difference/Legendre spectral method to solve the initial-boundary value time-fractional diffusion problem (1), on a finite domain, is considered. A finite element method for the time-fractional partial differential equation (1) on the sense of Riemann–Liouville is introduced in [10] and optimal order error estimates, both in semi-discrete and fully discrete cases, are obtained. Sidi Ammi and Jamiai have presented also a finite difference and Legendre spectral method for a time-fractional diffusion-convection equation for image restoration and a detailed error analysis was carried out [24]. In [25], Sidi Ammi and Torres consider a fractional nonlocal thermistor problem and develop a Galerkin spectral method. Some error estimates, in different contexts, are derived, showing that the combination of the backward differentiation in time and the Galerkin spectral method in space leads, for an enough smooth solution, to an approximation of exponential convergence in space [25]. Existence and uniqueness of solution for the fractional partial differential equation (1), with a left time Riemann–Liouville fractional derivative, is proved in [10] by using the Lax–Milgram Lemma. Here we propose a finite difference method in time and a finite element method in space to study the numerical solution of the time-fractional Caputo differential equation (1).

The outline of the paper is as follows. In Section 2, a finite difference scheme for solving the time-fractional diffusion equation is proposed, along with an unconditionally stability and convergence analysis. In Section 3, the finite element method is used and error estimates, in both time and space, are obtained. Then, some numerical tests are presented in Section 4, to verify the accuracy of the given method, comparing the obtained approximate results with the theoretical/exact ones. Some concluding remarks are given in Section 5. In the analysis of the numerical method that follows, we assume that problem (1) has a unique and enough regular solution.

2. Discretization in time: a finite difference scheme

In this section we consider the time discretization of (1). Define $A=-\Delta$ and $D(A)=H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$. Then the system (1) can be written, in abstract form, as

[TABLE]

Let $0=t_{0}<t_{1}<\cdots<t_{K}=T$ be a partition of $[0,T]$, where $t_{k}=k\Delta t$, $k=0,1,\ldots,K$, and $\Delta t=\frac{T}{K}$ is the time step. Following [16, 24], we discretize the Caputo derivative by a difference approach as follows: for all $0\leq k\leq K-1$,

[TABLE]

where $\tilde{R}_{k+1}$ is the truncation error satisfying

[TABLE]

and $c_{u}$ is a constant depending only on $u$. To continue the construction of the scheme, let us denote $b_{j}=(j+1)^{1-\alpha}-j^{1-\alpha}$, $j=0,1,\ldots,k$. It is easy to verify the following properties for $b_{j}$:

[TABLE]

Define the discretized fractional operator $L_{t}^{\alpha}$ by

[TABLE]

with

[TABLE]

Then,

[TABLE]

Let $t=t_{k+1}$. We can write (2) as

[TABLE]

Denote $u^{k}\approx u(t_{k})$ as the approximation of $u(t_{k})$. We define the following time stepping method:

[TABLE]

To complete the semi-discrete problem, we consider the boundary conditions

[TABLE]

and the initial condition

[TABLE]

We then obtain an equivalent form to (8):

[TABLE]

where $f^{k+1}=f(t_{k+1})$. For the particular case $k=0$, the scheme becomes

[TABLE]

If we define the error term $r^{k+1}$ by

[TABLE]

then it follows from (4) and (6) that

[TABLE]

Now we define some functional spaces endowed with standard norms and inner products that will be used in the remaining of the paper:

[TABLE]

[TABLE]

[TABLE]

where $L^{2}(\Omega)$ is the space of measurable functions whose square is Lebesgue integrable in $\Omega$. The inner products of $L^{2}(\Omega)$ and $H^{1}(\Omega)$ are defined, respectively, by

[TABLE]

while the corresponding norms are given by

[TABLE]

The variational weak formulation of equation (10) subject to the boundary condition (9) reads: find $u^{k+1}\in H^{1}_{0}(\Omega)$ such that

[TABLE]

$\forall v\in H^{1}_{0}(\Omega)$. Now we consider a stability result of the time discretization of equations (1).

Theorem 2.1**.**

Let $u^{k}$ be the approximation solution of (10). Then,

[TABLE]

Proof.

The result is proven by mathematical induction. First, when $k=0$, we have

[TABLE]

Then we get

[TABLE]

Note that $A$ is a positive definite elliptic operator, the eigenvalues of $A$ are $\lambda_{j}=j^{2}\pi^{2}$, $j=1,2,3\ldots$ It follows from the spectral method that the norm

[TABLE]

Hence, by using (15), we have

[TABLE]

Then,

[TABLE]

which suggests the result at the first step. Suppose now that the following hypothesis holds:

[TABLE]

We begin to prove that $\|u^{k+1}\|\leq\|u^{0}\|+\alpha_{0}(k+1)\|f\|_{L^{\infty}}$. From (10), we have

[TABLE]

Hence, by using (15) and (16), one has

[TABLE]

Finally, the last equality of (5) yields

[TABLE]

The proof is complete. ∎

We are now ready to prove error estimates in the $L^{2}$ norm for the error $u(t_{k})-u^{k}$ of the approximate solution $u^{k}$ of $u$, the exact weak solution $u(t_{k})$ of (8).

Theorem 2.2**.**

Let $u(t_{k})$ and $u^{k}$ be the solution of (7) and (8), respectively. Then,

[TABLE]

$k=1,2,\ldots,K$.

Proof.

We start by proving the following estimate:

[TABLE]

For that we use a standard induction procedure. Let $\varepsilon^{k}=u(t_{k})-u^{k}$. For $k=1$, we have, by calling together (7), (11) and (12), that the error equation is given by

[TABLE]

Hence, by using (15), we have

[TABLE]

With this in mind, and applying (13), one obtains that

[TABLE]

So, (17) is true for the case $k=1$. Suppose now that (17) holds for all $k=1,2,\ldots,K-1$. By gathering (7) and (10), we have

[TABLE]

It follows that

[TABLE]

Using the induction assumption, and the fact that the sequence $(b_{j})_{j}$ is decreasing, we obtain that

[TABLE]

Taking into account (5) in the above inequality, it follows that

[TABLE]

The auxiliary estimate (17) is then established. One can easily verify that $k^{-\alpha}b_{k-1}^{-1}\leq\frac{1}{1-\alpha}$ and $k\Delta t\leq T$, $k=1,2,\ldots,K$. Thus,

[TABLE]

and the proof is complete. ∎

In the coming section, we consider the space discretization of (1).

3. Discretization in space: a finite element scheme

The variational formulation of (1) consists to find $u(t)\in H_{0}^{1}(\Omega)$, such that

[TABLE]

More precisely, let $0=x_{0}<x_{1}<\cdots<x_{N}=1$ be an arbitrary space partition of $\Omega=[0,1]\subset\mathbb{R}$ and let $h=\max_{i}(x_{i+1}-x_{i})$. The set $\Omega$ can be a set of $\mathbb{R}^{n}$, $1\leq n\leq 3$. Let $S_{h}\subseteq H_{0}^{1}(\Omega)$ be a family of a finite element space consisting of piecewise linear continuous functions defined by

[TABLE]

Now consider the finite element method as follows: find $u_{h}(t)\in S_{h}$, such that

[TABLE]

Denote $A_{h}=-\Delta_{h}\leavevmode\nobreak\ :\leavevmode\nobreak\ S_{h}\leavevmode\nobreak\ \rightarrow\leavevmode\nobreak\ S_{h}$, which satisfies

[TABLE]

Let $P_{h}\leavevmode\nobreak\ :\leavevmode\nobreak\ H^{1}_{{}_{0}}(\Omega)\leavevmode\nobreak\ \rightarrow\leavevmode\nobreak\ S_{h}$ be the standard $L_{2}$ projection operator via the orthogonal relation

[TABLE]

and $R_{h}\leavevmode\nobreak\ :\leavevmode\nobreak\ H_{0}^{1}(\Omega)\leavevmode\nobreak\ \rightarrow\leavevmode\nobreak\ S_{h}$ be the elliptic or the Ritz projection defined by

[TABLE]

We can write (2) into abstract form as

[TABLE]

Denote by $u_{h}^{j}$ the approximation of $u(x,t_{j})$. We define the following time stepping method:

[TABLE]

Now consider the finite element discretization of problem (14) as follows: find $u^{k+1}_{h}\in S_{h}$, such that for all $v_{h}\in S_{h}$

[TABLE]

where the bilinear form $B_{h}(\cdot,\cdot)$ is defined by

[TABLE]

and the functional $F_{h}(v_{h})$ is given by

[TABLE]

Then, under enough regularity of the exact solution $u$, the following error estimate holds.

Theorem 3.1**.**

Let $u(t_{k})$ and $u_{h}^{k}$ be the solution of (2) and (19), respectively. Assume that $u\in H^{1}\left(H^{2}(\Omega)\cap H^{1}_{0}(\Omega),[0,T]\right)$. Then, the following inequality holds:

[TABLE]

Proof.

Let $\varepsilon_{h}^{k}=u_{h}^{k}-u(t_{k})$. We write

[TABLE]

where $\theta_{h}^{k}=u_{h}^{k}-R_{h}u(t_{k})$ and $\rho_{h}^{k}=R_{h}u(t_{k})-u(t_{k})$. We make use of the following inequality that can be find in [1, 4]:

[TABLE]

By the error estimate of the Ritz projection, we have

[TABLE]

In order to bound $\theta_{h}^{k}$, we use the error equation obtained from (19):

[TABLE]

From $\Delta_{h}R_{h}=P_{h}\Delta$ (see [27]), we hence obtain that

[TABLE]

where $\omega_{h}^{k+1}=\sigma_{h}^{k+1}+\tilde{R}_{k+1}$ and $\sigma_{h}^{k+1}=(I-R_{h})L^{\alpha}_{t}u(t_{k+1})$. Thus, we get

[TABLE]

Using the stability result of Theorem 2.1, we obtain that

[TABLE]

Hence, by using (22), we have

[TABLE]

Keeping in mind that $\|\tilde{R}_{k+1}\|\leq c_{u}\Delta t^{2-\alpha}$, we obtain that

[TABLE]

Together with this estimate, we get

[TABLE]

Hence, (21) is proved. ∎

4. Numerical validation

For completeness, our implementation is briefly described here.

4.1. Implementation

Considering problem (20), we express the function $u^{k+1}_{h}$ in terms of the finite piecewise linear elements, tent-line, global interpolation functions $\phi_{j}(x)$, $j=0,1,\ldots,N$,

[TABLE]

where $u_{j}(t_{k+1})=u(x_{j},t_{k+1})$ are unknowns of the numerical solution and $\phi_{j}$ are the global interpolation functions satisfying the cardinal interpolation property

[TABLE]

with $\delta_{ij}$ the Kronecker-delta symbol. By combining (20) and (23), and taking into account the homogeneous Dirichlet boundary condition $u^{k+1}_{0}=u^{k+1}_{N}=0$, we obtain the discrete system

[TABLE]

where

[TABLE]

[TABLE]

Since the matrix $\mathbb{M}+\alpha_{0}\mathbb{S}$ is symmetric positive definite, one can choose, for example, the conjugate gradient method to solve (24).

4.2. Numerical results

Now we present two numerical approximation examples to confirm our theoretical statements. The main purpose is to check the convergence behavior of the discrete solution with respect to the time step $\Delta t$ and the space step $\Delta x$ used in the computations.

Example 4.1**.**

Consider the time-fractional partial differential equation

[TABLE]

The right hand side $f$ and initial condition are selected as

[TABLE]

It is verified that the exact solution to the problem is

[TABLE]

The numerical results have been given by choosing $\Delta x=0.001$, $T=1$, and $\Delta t=0.01$, where $N=T/\Delta t$. Let $u^{n}$ denote the approximate solution, $u(t_{n})$ the exact solution, and $\varepsilon^{n}$ the error at $t=t_{n}$, that is, $\varepsilon^{n}=u^{n}-u(t_{n})$. Then we obtain Table 1 with the exact solution, the approximate solution, and the error for $\alpha=0.1,0.5,0.9$. We plot the exact solution, the approximate solution, and the error, for $\alpha=0.1,0.5,0.9$, in Figures 1, 2, and 3.

Example 4.2**.**

Consider the time-fractional partial differential equation

[TABLE]

with the forcing term and initial condition given by

[TABLE]

The exact solution is

[TABLE]

In this second example, we choose $\alpha=0.2$, $\Delta x=0.01$, $T=1$, $\Delta t=0.01$, and $N=T/\Delta t$. Let $u^{n}$ denote the approximate solution, $u(t_{n})$ the exact solution, and $\varepsilon^{n}=u^{n}-u(t_{n})$ the error at $t=t_{n}$. Figures 4(a) and 4(b) illustrate, respectively, the exact solution and the approximate solution at $t_{N}=1$. Figure 4(c) presents a plot of the error at $t_{N}=1$.

5. Conclusion

We have investigated a finite element method to Caputo time-fractional diffusion partial differential equations. A stability analysis is carried out and a convergent estimate is analyzed. We obtain error estimates in the $L_{2}$-norm between the exact solution and the approximate solutions in the fully discrete case. Two numerical examples are implemented and the numerical results are shown to be consistent with the theoretical results.

Acknowledgements

The authors were supported by the Center for Research and Development in Mathematics and Applications (CIDMA) of University of Aveiro, through Fundação para a Ciência e a Tecnologia (FCT), within project UID/MAT/04106/2019. They are very grateful to two anonymous referees, for careful reviews of their paper, and for the comments, corrections, and suggestions, which substantially helped them to improve the quality of the paper.