Approximate solutions of the $2$D space-time fractional diffusion equation via a gradient-descent iterative algorithm with Grünwald-Letnikov approximation

1.   Introduction

In modern days, fractional differential equations (FDEs) play an important role in applied mathematics, science, and engineering; see e.g. ^[1,2,3,4]. FDEs are powerful to model many physical phenomena due to a noninteger order in time and space.

Let us consider the standard two-dimensional ($2$D) diffusion equation (see e.g. ^[5]).

where $\mathbf{x} \in V \subseteq \mathbb{R}^2$ and $t \in [0, T] \subseteq \mathbb{R}^+$. Here, $\Delta$ denotes the Laplacian operator with respect to the spatial variable $\mathbf{x} = (x, y) \in \mathbb{R}^2$.

In physical contexts, the function $u(\mathbf{x}, t)$ represents the density of a diffusing material at location $\mathbf{x}$ and time $t$. Let us denote the Gamma function by $\Gamma$. From Eq (1.1), if we replace the first-order time derivative with the Caputo's derivative ^[2]

and replace the second-order space derivative with the Riemann-Liouville derivatives ^[6].

and similarly for $^{RL}D_y^{\alpha_2}u(x, y, t)$, then we get the $2$D space-time fractional diffusion equation (e.g. ^[7]).

Here, we impose the function $f(x, y, t)$ on the right hand side of the equation. Physical interpretations restrict the fractional orders so that $\alpha_1\in(1, 2]$, $\alpha_2\in(1, 2]$, and $\beta \in(0, 1]$. When $\beta = 1$, Eq (1.2) is known as the $2$D time fractional diffusion equation.

In recent years, there are many studies about theory of ($1$D or $2$D) space-time (or time) fractional diffusion equations and related equations, e.g. ^{[8,9,10,11,12,13]}. These differential equations were successfully used to describe many physical phenomena, e.g. random walk models ^{[14,15,16,17]}, and anomalous diffusion processes ^{[18,19,20,21]}.

In general, FDEs do not posses exact solutions in closed froms. Thus, numerical methods have been implemented for several types of FDEs, e.g. the variational iteration method ^[22], the homotopy analysis method ^[23,24], Adomian decomposition method ^[25,26], B-spline collocation schemes ^{[27,28,29,30]}, and the collocation method based on fractional Legendre functions ^[31]. There are also numerical methods using either finite difference or finite element methods to a discretization of certain FDEs; see e.g. ^{[32,33,34,35,36]}. Finite difference techniques based on Grünwald-Letnikov approximation for fractional derivatives were investigated by many authors, e.g. ^[1,37,38].

The main objective of the present paper is to propose an iterative algorithm to produce well-approximated solutions of the $2$D space-time fractional diffusion Eq (1.2), subjected to the zero Dirichlet boundary condition and the zero initial condition. First, we discretize Eq (1.2) via Grünwald-Letnikov approximation for the Riemann-Liouville and the Caputo derivatives (see Section 3). Then, we linearize the discretized equations into a linear system with a sparse cofficient matrix. In particular, we discuss the numerical scheme for its interesting special case, namely, the $2$D time fractional diffusion equation. To solve the linear system, we apply the gradient-descent technique to derive an iterative procedure with suitable directions and step sizes (see Section 4). We show that the produced approximate solutions converge to a unique solution for any given initial vector (see Section 5). Theoretical performance of the proposed algorithm are discussed through the convergence rate and error estimates. We verify the capability and theoretical performance of the proposed algorithm by making two numerical experiments (see Section 6). We compare the performance of the proposed algorithm with well-known iterative methods for linear systems in the literature, e.g. GI, LSI, SOR, and JOR algorithms.

2.   Preliminaries

In this section, we recall relevant background that used in later discussions.

2.1.   Grünwald-Letnikov approximation for fractional derivatives

The discretization of fractional derivatives is often done by finite difference schemes based on Grünwald-Letnikov type (see e.g. ^[37,38]). For the $2$D problem, we have the following approximation for Riemann-Liouville fractional derivatives of order $\alpha > 0$:

where the coefficients $g_{\alpha, k}$ are defined as

Alternatively, we have the recusive formula for $g_{\alpha, k}$:

We can compute the Grünwald-Letnikov approximation for Caputo's fractional derivatives via the difference formula.

Here, $\lfloor \cdot \rfloor$ denotes the floor function. When $\alpha \in (0, 1]$, the correction term on the right hand side of (2.3) is equal to zero when we consider the zero initial value, i.e., $u(x, y, t = 0) = 0$; see e.g. ^[38].

2.2.   Mittag-Leffler functions

The Mittag-Leffler function $E_{a, b}$ on two parameters $a, b > 0$ is a generalization of the exponential function defined by the Taylor series.

Mittag-Leffler functions arise naturally in certain fractional differential equations, in particular, in our numerical examples (Section 6).

2.3.   Auxiliary results from matrix analysis

All matrices and vectors considered throughout this paper are real. The Frobenius norm of a rectangular matrix $A$ is defined by $\Vert A \Vert_F = \sqrt{{\rm{tr}}(A^TA)}$. The Frobenius norm of a vector coincides with the Euclidean vector norm. The spectral norm $\Vert \cdot \Vert_2$ of a square matrix is defined to be the largest singular value of the matrix, or equivalently, the matrix norm induced by the Euclidean vector norm. Thus, we have the following relation (see e.g. ^[39]).

for any square matrix $A$ and vector $x$ of conformable sizes. The condition number $\kappa(A)$ of a square matrix $A$ is defined to be the ratio between the largest and the smallest singular values of $A$. In particular, if $A$ is invertible, we have

We recall gradient formulas for the trace of certain matrix products:

Lemma 2.1. (e.g. ^[39]) For matrices $A, B, X$ of compatible sizes, we have

3.   Discretization and linearization of the $2$D space-time fractional diffusion equation

In this section, we consider the $2$D space-time fractional diffusion equation, and its interesting special case, namely, the $2$D space fractional diffusion equation. We discretize these differential equations in which the Riemann-Liouville and the Caputo' derivatives are approximated by Grünwald-Letnikov approximation. Then we form a linear system from the discretized equations.

Let us consider the $2$D space-time fractional fractional diffusion equation of the form

on bounded domains $x\in [a_x, b_x]$, $y\in [a_y, b_y]$, and $t\in [0, T]$. The equation is subjected to the zero Dirichlet boundary condition

and the zero initial condition, i.e., $u(x, y, t = 0) = 0$.

3.1.   The $2$D space-time fractional diffusion equation

Consider Eq (3.1) with fractional orders $\beta \in (0, 1]$ and $\alpha_1, \alpha_2 \in (1, 2]$. We discretize $x_i = a_x+ih_x$, $y_j=a_y+jh_y$ and $t_k = kh_t$ where $h_x$, $h_y$ and $h_t$ are defined by

The approximated solution $u(x_i, y_j, t_k)$ at a point $(x_i, y_j, t_k)$ is written by $u_{ijk}$. For convenience, we define

From the difference Formula (2.3) and the Grünwald-Letnikov Formulae (2.1) and (2.2), we obtain

or equivalently,

for each $i\in \{1, \dots, N_x-1 \}$, $j\in \{1, \dots, N_y-1 \}$ and $k\in \{1, \dots, N_t\}$. For convenience, let us denote $N_{xy} = (N_x-1)(N_y-1)$ and $N = N_{xy} N_t$. Thus, Eq (1.2) can be linearized into a system of $N$ linear equations in $N$ unknowns $u_{111}, \dots, u_{N_x-1, N_y-1, N_t}$.

To form a linear system, let

Then put $\mathbf{u} = {\rm{Vec}}(U)$ and $\mathbf{f} = \sigma {\rm{Vec}}(F)$. Here, ${\rm{Vec}}(\cdot)$ is an operator that turns a matrix into a column vector by stacking columns of the matrix consecutively; see e.g. [40, Ch. 4]. Hence, we obtain a linear system

where the coefficient $P$ is a block lower-triangular matrix

and

We shall seek for a well-approximated solution $\mathbf{u}$ of the linear system (3.3). Note that the system (3.3) has a unique solution if and only if $P$ is invertible, or equivalently, $A$ is invertible. Once we find the vector $\mathbf{u}$, we can obtain the matrix $U$ due to the injectivity of the operator ${\rm{Vec}}(\cdot)$.

3.2.   An interesting special case: The 2D space fractional diffusion equation

Here, we consider the $2$D space fractional diffusion equation

where $\alpha_1, \alpha_2 \in (1, 2]$, the spatial domain is $[a_x, b_x]\times [a_y, b_y]$, and the temporal domain is $[0, T]$. Eq (3.4) is a special case of Eq (3.1) for which the fractional order $\beta$ is equal to $1$.

We discretize Eq (3.4) by computing an approximated solution at $(x_i, y_j, t_k)$ with $x_i = a_x+ih_x$, $y_j = a_y+jh_y$ and $t_k = kh_t$ where $h_x$, $h_y$ and $h_t$ are defined as (3.2). The first-order and the fractional-order derivatives are approximated by the forward time difference method and the Grünwald-Letnikov approximation (2.1) and (2.2), respectively. For convenience, we define

From Eq (3.4), we obtain that for each $i\in \{1, \dots, N_x-1 \}$, $j\in \{1, \dots, N_y-1 \}$ and $k\in \{1, \dots, N_t\}$,

or equivalently,

We use the notations $N_{xy}$, $N$, $U$, $F$, $\mathbf{u}$, and $\mathbf{f}$ as same as those in Subsection 3.1. Then we can put the discretized equations into a system of $N$ linear equations in $N$ variables $u_{111}, \dots, u_{N_x-1, N_y-1, N_t}$ and displayed as

In this case, $P$ is a block lower-triangular matrix as follows:

where

4.   Proposing a gradient-descent iterative algorithm

In this section, we apply the gradient-descent technique to derive an iterative procedure for solving the linear system (3.3) and, in particular, the system (3.5). Assume that the coefficient matrix $P$ is invertible. Then we can solve for the analytical solution $\mathbf{u}^*$ directly by

For approximate solutions, we measure an error using the quadratic norm-error function

We introduce a gradient-descent iterative solver described by the recursive formula

We pick an arbitrary initial solution $\mathbf{u}_0$, so that Eq (4.3) iteratively computes the next solution $\mathbf{u}_{k+1}$ and form a sequence of approximate solutions. Due to the gradient direction $\nabla \Omega(\mathbf{u}_k)$ and a convergence factor $\tau_{k+1}$, the sequence $\{ \mathbf{u}_k\}_{k = 0}^{\infty}$ would converge to the exact solution $\mathbf{u}^*$.

To find the gradient of $\Omega$, we apply Lemma 2.1 to derive

Hence, the iterative equation takes the form

According to the gradient-descent, the convergence factor $\tau_{k+1}$ is generated in order to minimize the error occurring at each iteration. Thus, we define a new function $s_{k+1} : [0, \infty)\to \mathbb{R}$ to be an error at step $k+1$:

Using properties of the matrix trace, we obtain the derivative of $s_{k+1}$ as follows:

It is easy to check that the second-order derivative $\frac{d^2}{d\tau^2}s_{k+1}(\tau)$ is positive. Setting $\frac{d}{d\tau}s_{k+1}(\tau) = 0$, we have the minimizer of the function $s_{k+1}(\tau)$ as follows:

We call $\{ \tau_{k+1}\}_{k = 0}^{\infty}$ the sequence of optimal convergence factors. We now describe the Frobenius norms $\Vert P^T(\mathbf{f}-P \mathbf{u}_k) \Vert_F$ and $\Vert PP^T(\mathbf{f}-P \mathbf{u}_k) \Vert_F$ more precisely. To avoid duplicated multiplications, let $\mathbf{y} = P^T \mathbf{f}$ and $V = P^T P$. Consider

Let us denote $\mathbf{y}_i$ and $v_{ij}$ for the $i$th entry of $\mathbf{y}$ and the $(i, j)$th entry of $V$, respectively. It follows that

In a similar way, by denoting $\mathbf{z} = P\mathbf{y}$ and $W = PV$, we obtain

We combine the direction and the step size altogether to obtain the following iterative algorithm for solving Eqs (3.3) and (3.5), see Algorithm 1.

To break the procedure, one can impose a stopping rule $\Vert P \mathbf{u}_k- \mathbf{f} \Vert_F < \epsilon$ or $\Vert \mathbf{u}_{k}- \mathbf{u}_{k-1} \Vert_F < \epsilon$, where $\epsilon$ is a satisfactory error. Note that the coefficient matrix $P$ is sparse, so that the computational procedures need not so much time.

5.   Capability and theoretical performance of the proposed algorithm

In this section, we investigate the capability and theoretical performance of Algorithm 1 through error esimates and the convergence rate.

To show that Algorithm 1 is applicable for any initial vector, let us recall the following approximations for strongly convex functions on $\mathbb{R}^n$. Recall also that for real symmetric matrices $A$ and $B$ of the same size, the matrix ordering $A \leqslant B$ means that the difference $B-A$ is positive semidefinite.

Lemma 5.1. (^[41]) Let $f: \mathbb{R}^n\to \mathbb{R}$ be a strongly convex function, i.e., there exist two numbers $\phi, \Psi \geqslant 0$ such that the matrix orderings $\phi I \leqslant \nabla^2f(x) \leqslant \Psi I$ hold for all $x \in \mathbb{R}^n$. Then, for any $u, v \in \mathbb{R}^n$, we have

Theorem 5.2. Consider the linear system $(3.3)$ where the coefficient matrix $P$ is invertible. Denote by $\kappa$ the condition number of $P$.Then the following statements hold:

(i) For any initial vector $\mathbf{u}_0$, the sequence of approximate solutions $\{ \mathbf{u}_k\}$ generated by Algorithm 1 converges to a unique solution $\mathbf{u}^*$.

(ii) Error estimates of $\Vert P \mathbf{u}_k-\mathbf{f} \Vert_F$ compared to the preceding step and the initial step are described by the following inequalities:

In particular, the relative error decreases at every iteration.

(iii) Algorithm 1 has the convergence rate (in regard to the relative error $\Vert P \mathbf{u}_k-\mathbf{f} \Vert_F$) governed by $\sqrt{1-\kappa^{-2}}$.

(iv) Moreover, the following error estimates for $\Vert \mathbf{u}_k- \mathbf{u}^* \Vert_F$ hold:

Proof. If there exists an integer $k > 0$ so that a vector $\mathbf{u}_k$ makes $\nabla \Omega(\mathbf{u}_k)$ to be zero, we have $\mathbf{u}_k = \mathbf{u}^*$ and the result holds. So we assume that $\nabla \Omega(\mathbf{u}_k)\neq 0$ for all $k$. Now, we compute the second-order derivative $\nabla^2 \Omega(\mathbf{u}) = P^T P$, which is a positive semidefinite matrix. Let $\lambda_{\min}$ and $\lambda_{\max}$ be the smallest and the largest eigenvalues of $P^T P$, respectively. From the spectral theory of matrices, we obtain

where $I$ is the identity matrix of compatible size. Hence, $\Omega$ is a strongly convex function in which $\phi = \lambda_{\min}$ and $\Psi = \lambda_{\max}$. Applying Eq (5.2) to the function $s_{k+1}(\tau)$, we have

Let us define the right-hand side (RHS) of (5.7) to be a function of $\tau$, namely,

We can verify that $1/ \lambda_{\max}$ is a minimizer of $f$. Minimizing (5.7) both sides by $\tau = 1/ \lambda_{\max}$, we obtain

Applying Eq (5.1) to the function $\Omega {k+1}(\tau)$, we have

Similarly, we have that $\tau = 1/ \lambda_{\min}$ is a minimizer of the RHS of (5.9). We continue in this fashion to obtain

Now, $\Vert\nabla \Omega (\mathbf{u}_k) \Vert^2_F\geqslant 2 \lambda_{\min} \Omega (\mathbf{u}_k)$, and hence by taking account of (5.8), we get

Since $P$ is an invertible matrix, all eigenvalues of $P^T P$ are strictly positive, and thus $1-\kappa^{-2} < 1$. By recurring the above inequality, it follows that for any $k \in \mathbb{N}$,

Therefore, $\Omega (\mathbf{u}_k) \to 0$, i.e., $\mathbf{u}_k \to \mathbf{u}^*$ as $k\to\infty$. Moreover, from the error bounds (5.10) and (5.11), we get the inequalities (5.3) and (5.4), respectively. Thus, the asymptotic behavior of $\mathbf{u}_k$ with respect to the relative error $\Vert P \mathbf{u}_k- \mathbf{f} \Vert_F$ is governed by $\sqrt{1-\kappa^{-2}}$.

From the estimate (5.4) and the norm properties (2.4) and (2.5), we derive

Similarly, we can have the error estimate $\Vert \mathbf{u}_k- \mathbf{u}^* \Vert_F$ compared to the preceding step $\Vert \mathbf{u}_{k-1}- \mathbf{u}^* \Vert_F$ using Eqs (2.4) and (2.5) together with the bound (5.3).

Theorem 5.2 implies that the condition number $\kappa$ of the sparse matrix $P$ effects the convergence behaviour. Indeed, Algorithm 1 converges faster to the exact solution if $\kappa$ is close to $1$.

6.   Capability and performance of the algorithm through numerical experiments

In this section, we implement numerical experiments to perform the capability and performance of Algorithm 1, which is denoted by TauOpt. All experiments have been carried out by MATLAB R2020b with PC environment Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz, RAM 8.00 GB. To measure the time consuming (in seconds) taken for each simulation, we apply the tic and toc functions in MATLAB and abbreviate them TC. At the $k$th iteration step, we concern the relative error $\Vert \mathbf{f}-P \mathbf{u}_k \Vert_F$. In each example, we provide tables and figures to demonstrate the accuracy and efficiency of algorithms. We compare the performance of the proposed algorithm with famous iterative methods for the linear system, for instance, GI algorithm ^[42], LSI algorithm ^[42], SOR algorithm, and JOR algorithm ^[43].

Example 6.1. Consider the $2$D space-time fractional diffusion equation

where

The domains are $x\in[0, \pi]$, $y\in[0, \pi]$, and $t\in[0, 0.1]$. Indeed, the exact solution for this equation subjected to the zero Dirichlet boundary condition and the zero initial condition is given by

According to the numerical scheme explained in Subsect. 3.1, we discretize Eq $(6.1)$ by choosing partition numbers $N_x = 10$, $N_y = 10$, and $N_t = 20$. Then the coefficient matrix $P$ of the associated linear system $(3.3)$ is of size $1620\times 1620$. In this case, the condition number of $P$ is $\kappa = 10.2976$. Theoretically, Algorithm 1 has the convergence rate governed by $\sqrt{1-\kappa^{-2}} = 0.9953$. Table 1 shows the comparisons between numerical solutions and the exact solution for certain values of $x, y, t$ in the domain. After running $250$ iterations, we see that the $4$-digits numerical solutions are very close to the exact solution. Also, Figure 1 illustrates the $3$D plot of the exact and numerical solutions in the cases that $t = 0.01, 0.02, \dots, 0.06$. Besides, we compare the performance of Algorithm 1 with SOR, JOR, GI and LSI algorithms. The results, shown in Table 2 and Figure 2, indicate that Algorithm 1 performs well in both iteration numbers and computational time.

Example 6.2. Consider the $2$D space fractional diffusion equation

where

subjected to the zero Dirichlet boundary condition and the zero initial condition. The domains are $x\in[0, \pi]$, $y\in[0, 1]$, and $t\in[0, 0.1]$. Indeed, the exact solution for this problem is given by

According to the discussions in Subsection 3.2, we discretize Eq $(6.2)$ using partition numbers $N_x = 10$, $N_y = 10$ and $N_t = 20$. In this case, the condition number of the coefficient matrix $P \in \mathbb{R}^{1690\times 1690}$ is $\kappa = 26.7484$. Hence, Algorithm 1 has the convergence rate governed by $\sqrt{1-\kappa^{-2}} = 0.9993$. We run Algorithm 1 for $500$ iterations and report the comparison between numerical solutions and the exact solution in Table 3. Figure 3 illustrates the $3$D plot of the exact and numerical solutions in the cases that $t = 0.01, 0.02, \dots, 0.06$. The results from Table 3 reveal that numerical solutions are very close to the exact solution.

In addition, we compare the performance of Algorithm 1 (TauOpt) with SOR, JOR, GI, and LSI algorithms. The results of running $500$ iterations are shown in Table 4. Figure 4 displays the logarithmic relative error of each algorithm. It is seen that Algorithm 1 performs well in both iteration numbers and time consuming.

Note that in Examples 6.1 and 6.2, the condition numbers $\kappa$ seem to be too much, so that the convergences of the proposed algorithm would be not fast. However, due to the fact that the coefficient matrix $P$ is sparse and the convergence factor is chosen by an optimazation technique, the iterative procedures lead to a desire solution with a satisfactory error in a short time.

7.   Conclusions

We discretize the $2$D space-time fractional and the $2$D space fractional diffusion equations via the finite difference scheme of Grünwald-Letnikov approximation. We transform the discretized equations into a sparse linear system with coefficient matrix $P$. We propose an iterative solver based on the technique of gradient-descent. The analysis confirms the capability of the proposed algorithm as long as the matrix $P$ is invertible with the convergence rate governed by $\sqrt{1-\kappa^{-2}}$, where $\kappa$ is the condition number of $P$. The numerical experiments indicate the accuracy and the efficiency of the proposed algorithm compared to well-known iterative solvers for linear systems. Indeed, the iterative procedures of the proposed algorithm lead to a desire solution with a satisfactory error in a short time since the coefficient matrix $P$ is sparse, the convergence factor is chosen by an optimazation technique, and the procedures avoid duplicated computations.

Conflict of interest

All authors declare no conflicts of interest in this paper.