High order compact difference scheme for solving the time multi-term fractional sub-diffusion equations

1.   Introduction

Fractional differential equations appear in many practical and scientific applications, such as finance, engineering, chemistry, physics ^{[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,39]}. In general, it is difficult to find analytic solution of most fractional differential equations, thus, it is necessary to use numerical techniques to solve fractional differential equations.

Generally, Riemann-Liouville fractional derivative and the Caputo fractional derivative are most favorable definitions of fractional derivative. Based on the interpolation approximation, the commonly used approximation formulas is $L1$ formula, which is derived from a piecewise linear interpolation approximation. In ^[25], the numerical accuracy of this formula is proved to be $(2-\alpha)$ $(0 < \alpha\leq1)$ order. Using the piecewise quadratic interpolation approximation, the $L2-1$ formula for Caputo fractional derivative ${_{\; 0}^{C}}{\cal D}_{t}^{\alpha}f(t)$ is derived in Gao et al.^[26]. In order to obtain the higher order of the Caputo fractional derivative, Alikhanov ^[27] derived a $L2-1_{\delta}$ formula to approximate Caputo fractional derivative at a superconvergent point. In the above works, the proposed method is only devoted to the numerical approximation of the one-term fractional derivative. In fact, the multi-term fractional differential equations has been proved to be valuable models for describing many processes in practice ^{[28,29,30,31,32]}.

In order to find accurate and efficient numerical algorithms for multi-term fractional differential equations, some research work on this subject has been made. For example, the collocation method in ^[30], the iterative method ^[29], the fractional predictor-corrector method in ^[33] etc. In ^[34], the multi-term variable-distributed order diffusion equation is studied by $L1$ formula for the approximation of time-fractional derivatives. Until now, there is relatively little discussion on high order numerical method for the time multi-term fractional sub-diffusion equations.

In this paper, we consider the following time multi-term fractional sub-diffusion equations:

where $\rho_0, \rho_1, \ldots, \rho_m$ are some positive constants, $0\leq \alpha_m < \alpha_{m-1} < \ldots < \alpha_0\leq1$ and the term ${_{\; 0}^{C}}{\cal D}_{t}^{\alpha}u(x, t)$ represents the Caputo fractional derivative of order $\alpha$, which is defined by

2.   Notations and preliminary lemmas

Let $h = L/M$ be the spatial step, $\tau = T/N$ be the time step, which $M, N$ are positive integers. We partition $[0, L]$ into a mesh by the mesh points $x_{i} = ih$ $(0\le i\le M)$. Denote $t_{n} = n\tau$ $(0\leq n\leq N)$. Let $u(x, t)$ be the solution of the problem (1.1). Define the grid functions

Define the compact operator^[36]

and

where $k'$ and $k''$ denote the first and second order derivative of the coefficient function $k(x)$, respectively.

For the requirement of analysis, we assume the there exist positive constants $c_1$ and $c_2$ such that

Let ${\cal S}_h = \{u = (u_0, u_1, \cdots, u_M), u_0 = u_M = 0\}$ be the space of the grid functions defined in the spatial mesh. For any grid functions $u, v\in{\cal S}_h$, we define the inner product and norm as follows:

For any grid functions $u, v\in{\cal S}_h$, according some simple calculations, we have^[38]

By the following Lemma 2.5, we can see that $\left ({\cal H}u, u\right)\ge \frac{1}{4}\|u\|^2$, when the spatial size $h$ is sufficient small. For convenience, we introduce the discrete inner product and the norm as follows:

As done in ^[35], we denote

Lemma 2.1. ^[35]When $m\ge1$, the Newton iteration sequence $\{\sigma_p\}_{p = 0}^{\infty}$, generated by

is monotonically decreasing and convergence to $\sigma^*$ which is a unique positive root of $G(\sigma) = 0$ (Lemma 2.1 in ^[35]).

For $0 < \alpha < 1$, we define

Let

when $n\ge 2$, denote^[38]

Denote $t_{n-1+\sigma} = (n-1+\sigma)\tau$ and consider $\sum\limits_{r = 0}^m\rho_r{_{\; 0}^{C}}{\cal D}_{t}^{\alpha_r}u(x_i, t_{n-1+\sigma})$ for any $\alpha\in[0, 1]$. The following lemma gives a numerical differentiation formula to approximate $\sum\limits_{r = 0}^m\rho_r{_{\; 0}^{C}}{\cal D}_{t}^{\alpha_r}(x, t)$ at the point $t = t_{n-1+\sigma}$ and reveals its numerical accuracy.

Lemma 2.2. Suppose $u(x, t)\in C^{(4, 3)}([0, L]\times[0, T])$, it holds that

where

Proof. It follows from Theorem 2.1 of ^[35].

Lemma 2.3. ^[35] Given any non-negative integer $m$ and positive constants $\rho_0, \rho_1, \cdots, \rho_m$, for any $\alpha_i\in[0, 1], i = 0, 1, \cdots, m$, where at least one of $\alpha_i$'s belongs to $(0, 1)$, it holds

Applying the Lemmas 2.2 and 2.3 in ^[36] and following the similar procedure in ^[37], we have the following lemma.

Lemma 2.4. It holds that

and

where $C^*$ is a positive constant.

Proof. The proof follows from (31)–(33) in ^[37].

Lemma 2.5. For any grid function $u\in{\cal S}_h$, we have

Proof. Using the definition of compact operator ${\cal H}$, we obtain

Taking the inner product of the above equation with $u^n$ and applying the discrete Green formula, we obtain

Noticing (2.2), we have $|u|_1^2\le\frac{4}{h^2}\|u^n\|^2$, thus, we get

According to some simple process, we have

Noting (2.1) and $u_M^n = 0$, we obtain

This proves (2.3).

Lemma 2.6. Assume $u = \{u_{i}^{n}\; |\; 0\le i\le M, \; 0\leq n\leq N\}\in {\cal S}_{h}$, we obtain

Proof. By the discrete inner product $\langle \cdot, \cdot\rangle$, we have

Using the Corollary 1 of ^[27], we have

The (2.4) follows immediately.

3.   Compact difference schemes for the time multi-term fractional sub-diffusion equation

Based on the above lemmas, we now discrete (1.1) into a compact finite difference system. Considering the governing equation of (1.1) at the point $(x_i, t_{n-1+\sigma})$, we have

It follows from Lemmas 2.2 that

where

Using Lemma 2.3, we have

where

Then we obtain the compact difference scheme as follows:

Theorem 3.1. Assume $u(x, t)$ is the solution of problem (1.1). The truncation error $(R_{tx})_i^n$ satisfies

where $\mathcal{C}$ is a positive constant.

Theorem 3.2. The difference scheme (3.1) is uniquely solvable.

Proof. It is clear that the value of $u^0$ is determined. Assume $\{u^k|0\le k\le n-1\}$ has been determined, then we obtain a linear system of equations with respect to $u^n$. Considering the corresponding homogeneous system:

Multiplying the governing equation of (3.2) by $u_i^n$ and summing up for $i$ from $1$ to $M-1$, we obtain

Applying Lemma 2.4, we have

Let $\varepsilon = \frac{2}{3}-\frac{h^2}{24}c_2$, when the spatial step $h$ is sufficient small such that $\varepsilon > 0$, thus, we get

Noticing (3.2), we can get $u^n = 0$. This completes the proof.

We now carry out a prior estimate of the difference scheme (3.1).

Theorem 3.3. Suppose $\{u_i^n|0\le i\le M, 1\le n\le N\}$ is the solution of the difference scheme (3.1). Then when $\tau^{\alpha_0}\leq\frac{2c_{0, n}^{(\alpha)}}{\Gamma(2-\alpha_0)\sigma^2}$, we have

Proof. At first, we take an inner product of (3.1) with $\sigma u^n+(1-\sigma)u^{n-1}$ yield

By Lemma 2.5, we have

Substituting (3.5) into (3.4), we have

Using the Cauchy-Schwarz inequality, we obtain

Substituting (3.7) into (3.6), we can get

or equivalently,

Since $\|u^n\|_{\epsilon}^2\ge\frac{1}{4}\|u^n\|^2$ and $\sigma\in[\frac{1}{2}, 1)$(Lemma 2.1 in ^[35]), we have

Thus we obtain

When $\tau^{\alpha_0}\leq\frac{2c_{0, n}^{(\alpha)}}{\Gamma(2-\alpha_0)\sigma^2}$, we have

Denote

It follows from the above inequality that

We prove by induction that

This completes the proof of (3.3).

4.   Applications and numerical results

In this section, we are devoted to some numerical illustration on the theoretical results in the previous sections. Denote the norm errors as follows:

The temporal and spatial convergence orders are computed, respectively, by

Example 4.1. We first consider a problem with $\beta_{1}\not = 0$. This problem is governed by the Eq (1.1) in $[0, 1]\times [0, 1]$ with $k(x, t) = \frac{1}{1+x+t}$, $q(x, t) = (x+2t)^2$ and

The boundary and initial conditions are given with

It is easy to check that $v(x, t) = e^{t-2x}$ is the solution of this problem.

At first, we test the temporal convergence order of the compact difference scheme (3.6) for different $\alpha$. In this test, we let the spatial step $h = 1/1000$. Tables 1 and 2 give respectively the errors ${\rm E}_{\nu}(\tau, h)$ $(\nu = 1, 2, \infty)$ and the temporal convergence orders ${\rm O}_{\nu}^{\rm t}(\tau, h)$ $(\nu = 1, 2, \infty)$. As expected from our theoretical analysis, the computed solution $u_{i}^{n}$ has the second-order temporal accuracy.

In Figures 1 and 2, we plot the errors ${\rm E}_{\nu}(\tau, h)$ $(\nu = 1, 2, \infty)$ as a function of the time step $\tau$ with the spatial step $h = 1/1000$ for $\alpha = (1/3, 1/4, 1/5)$ and $\alpha = (1, 1/2, 0)$. In these two figures, logarithmic coordinates are considered. We see from Figures 1 and 2 that the variation of the error with the time step $\tau$ is linear and parallel to the line $y = 2x$. This numerically confirms the temporal $2$th-order convergence of the numerical solution $u_{i}^{n}$.

We next compute the spatial convergence order of the compact difference scheme (3.6). Tables 3 and 4 present the errors ${\rm E}_{\nu}(\tau, h)$ $(\nu = 1, 2, \infty)$ and the spatial convergence orders ${\rm O}_{\nu}^{\rm s}(\tau, h)$ $(\nu = 1, 2, \infty)$ with the time step $\tau = 1/20000$. Tables demonstrate that the compact difference scheme (3.6) has the fourth-order spatial accuracy.

Figures 3 and 4 show that the variation of the error with the spatial step $h$ is parallel to the line $y = 4x$. This corresponds to the spatial fourth-order convergence of the numerical solution $u_{i}^{n}$.

5.   Conclusions

We have presented and analyzed a high-order compact finite difference method for a class of time multi-term fractional sub-diffusion equations. In our method, a higher accurate interpolation approximation for a linear combination of the multi-term fractional derivatives in the Caputo sense and a fourth-order compact finite difference approximation is used for the spatial derivative. We have proved that the resulting scheme is unconditionally stable and convergent. We have also provided the optimal error estimates in the discrete $H^{1}$, $L^{2}$ and $L^{\infty}$ norms. The error estimates show that the proposed method has the second-order temporal accuracy and the fourth-order spatial accuracy. Numerical results confirm our analysis and demonstrate the efficiency of our method.

Acknowledgments

The authors would like to thank the referees for their valuable comments and suggestions which improved the presentation of the paper. This work was supported by National Natural Science Foundation of China (No.11401363), by the Key Scientific Research Projects of Colleges and Universities in Henan Province (No.21A110020).

Conflict of interest

This work does not have any conflict of interest.