Effective difference methods for solving the variable coefficient fourth-order fractional sub-diffusion equations

1.   Introduction

In this paper we consider the numerical approximations of the following problem

where ${\kappa}\geq 0$ is given constant, $\varphi (x), {\alpha_1}(t), {\alpha_2}(t), {{\beta}_1}(t), {{\beta}_2}(t)$ and $f(x, t)$ are given sufficiently smooth functions satisfying $\varphi (0) = {{\alpha}_1}(0), \varphi (L) = {{\alpha}_2}(0), \varphi ''(0) = {{\beta} _1}(0)$ and $\varphi ''(L) = {{\beta}_2}(0)$, $\; _0^CD_t^{\alpha}u(x, t)$ denotes Caputo fractional derivative defined by

And we suppose that there exist two constants $C_1$ and $C_2$ such that $0 < {C_1} \le {\omega}(x) \le C_2$ for $0\le x \le L$.

More and more attention has been paid to the fractional differential equations (FDEs) due to its application foreground in chemistry, physics, finance and hydrology in the past twenty years ^[1,2,3,4]. As we know, the analytic solutions of FDEs are very difficult to obtain, some efficient numerical methods should be considered, especially fast algorithms with high order accuracy. Some essential definitions and properties of fractional derivatives can refer to monograph ^[5].

This target problem in Eq (1.1) is frequently employed to simulate some phenomena in physics, such as wave propagation in beams, brain warping, ice formation and designing special curves on surfaces and so on, e.g., ^{[6,7,8,9,10,11]} and their references.

Up to now considerable works have been done from theoretical and numerical point of view for fourth-order fractional diffusion equations. For instance, Hu and Zhang successively presented a finite difference scheme for the fourth-order fractional diffusion-wave and sub-diffusion equations, and a compact difference scheme for the former, see ^[12,13]. Ji et al.^[14] constructed a compact difference scheme for the fourth-order fractional sub-diffusion equation under the fist Dirichlet boundary conditions. Zhang and Pu ^[15] presented a compact difference scheme for such equation by $\mathcal{L}2 - 1_{\sigma}$ formula ^[16]. Ran and Zhang ^[17] presented a new compact difference schemes for the such equation of the distributed order.

However, most of the work focus on the constant coefficient case. Recently, Zhao and Xu ^[18] presented a compact difference scheme for the time fractional sub-diffusion equation with the variable coefficient under the Dirichlet boundary conditions. Subsequently, based on the subtle decomposition of the coefficient matrices, Vong, Lyu and Wang ^[19] presented a compact difference scheme to solve the equations under Neumann boundary conditions. But the above works has only accuracy of order $2-\alpha$ in time.

In this paper, our attention will be paid on the higher order difference scheme for solving the variable coefficient equations under the second Dirichlet boundary conditions For this purpose, we use the $\mathcal{L}2 - 1_{\sigma}$ formula to approximate the Caputo fractional derivative. Unlike the integer order case, the time fractional derivative requires all history information. In order to reduce the computational complexity, we also construct a fast difference scheme. The stability and convergence of both schemes are proved in detail.

The structure of this paper is as follows: In Section 2, some necessary notations and lemmas are first introduced and a second-order difference scheme for the target problem (1.1)–(1.4) is constructed. In Section 3, an important priori estimate is first proved, and the unconditional stability and convergence of scheme are obtained. In Section 4, a fast second-order difference scheme is presented, and the corresponding unconditional stability and convergence are also strictly proved.

In Section 5, a difference scheme based on nonuniform time grids is first presented, and some numerical examples are provided to verify the theoretical results. A brief conclusion is given finally.

2.   Derivation of the $\mathcal{L}2 - 1_{\sigma}$ scheme

Let $h = L/M$ and $\tau = T/N$, where M, N are two positive integers. Denote $x_i = ih, 0\le i\le M, t_n = n\tau, 0\le n\le N$, $\Omega_h = \{x_i\mid0\le i\le M\}, \; \Omega_{\tau} = \{t_n\mid0\le n\le N\}$. Let $\mathcal{V}_h = \{v\mid\; v = (v_0, v_1, \cdots, v_M)\}$ be grid function space on $\Omega_h$, and $\mathcal{\mathring{V}}_h = \{v\mid\; v\in\mathcal{V}_h, v_0 = v_M = 0\}$. Also we denote $\sigma = 1-\frac{\alpha}{2}, t_{n+{\sigma}} = (n+{\sigma})\tau$ and $\omega(x_i) = \omega_i$.

For $u\in \mathcal{V}_h$, we define

For any $u, v\in\mathcal{\mathring{V}}_h$, we define the inner products

and norms

In ^[16], Alikhanov developed a new second order difference formula (called $\mathcal{L}2 - {1_{\sigma}}$ formula) for the Caputo fractional derivative, which can be expressed in the following lemma.

Lemma 2.1 (^[16]). Suppose $\alpha\in (0, 1), {\sigma} = 1-\frac{{\alpha}}{2}$ and $u(t)\in {C^3}[0, T]$. It holds

where

in which $C_0^{(n)} = {a_0} = {{\sigma}^{1-{\alpha}}}$ for $n = 1$, and

for $n\ge 2$, where ${a_j} = {(j+{\sigma})^{1-{\alpha}}}-{(j-1+{\sigma})^{1-{\alpha}}}$ and ${b_j} = \frac{1}{{2-{\alpha}}}[{(j+{\sigma})^{2-{\alpha}}}-{(j-1+{\sigma})^{2-{\alpha}}}]-\frac{1}{2}[{(j+{\sigma})^{1-{\alpha}}}+{(j-1+{\sigma})^{1-{\alpha}}}]$ for all $j \ge 1$.

Let $v(x, t) = \frac{{{{\partial}^2}u}}{{{\partial}{x^2}}}$. Then the problem Eqs (1.1)–(1.4) can be written in the equivalent system

Suppose $u(x, t)\in C_{x, t}^{(6, 3)}([0, L]\times[0, T])$. Define

Considering the Eqs.(2.1)–(2.2) at the point $({x_i}, {t_{n-1+{\sigma}}})$, we obtain

Using Taylor expansion

where $U_i^{n-1+{\sigma}} = {\sigma}U_i^n+(1-{\sigma})U_i^{n-1}$. Then we obtain

and

Using Lemma 2.1, it follows from Eq (2.5), Eq (2.6) that

and there exists a constant $C_r$ such that

Omitting the small terms $({R_1})_i^n$ and $({R_2})_i^n$ in Eq (2.7) and Eq (2.8), we present the difference scheme (called $\mathcal{L}2-{1_\sigma}$ scheme) for the equivalent system (2.1)–(2.4) as follows

where the initial-boundary conditions Eq (2.3), Eq (2.4) have been used.

Theorem 2.2. The above difference scheme (2.10)–(2.13) is equivalent to

Proof. Since

It follows from Eq (2.11) and Eq (2.13) that

This together with Eq (2.10), we get Eq (2.14) and Eq (2.16). Eq (2.15) can be obtained by substituting Eq (2.11) into Eq (2.10). This proof is completed.

The above equivalent form Eqs (2.14)–(2.18) will be used only in calculation.

3.   Solvability, stability and convergence of the $\mathcal{L}2-{1_\sigma}$ scheme

We first introduce the following essential lemmas.

Lemma 3.1 (^[16]). Suppose ${\alpha} \in (0, 1)$ and $C_k^{(n)}$ is defined in Lemma 2.1. It holds that

Lemma 3.2 (^[16]). Suppose $u = \{ {u^n}\mid0 \le n \le N\}$ is a grid function defined on ${{\Omega}_\tau }$. It holds that

Lemma 3.3 (^[20,21]). For any $u \in \mathcal{\mathring{V}}_h$, it holds that

The following Lemma will be used in the analysis of the difference scheme.

Lemma 3.4. For any $u \in \mathcal{\mathring{V}}_h$, it holds that

Proof. The proof is straightforward from the definition of $||\cdot||$ and $||\cdot||_{\omega}$.

We next show the priori estimate of the scheme (2.10)–(2.13).

Theorem 3.5. Suppose $\{w_i^n\mid0\le i\le M, 0\le n\le N\}$ and $\{z_i^n\mid0\le i\le M, 0\le n\le N\}$ satisfy the following difference scheme

Then, it holds that

Proof. Taking the inner product of Eq (3.1) by $w^{n-1+{\sigma}}$, we get

Taking the inner product of Eq (3.2) by ${{\omega}{z^{n-1+{\sigma}}}}$, we get

From Eq (3.6) and Eq (3.7), it yields that

Applying the discrete Green formula gives that

Substituting Eq (3.9) into Eq (3.8), we obtain

From Eq (3.2), we have

Multiplying Eq (3.11) by $h\omega_i$ and summing up for $i$ from $1$ to $M-1$, we get

Substituting Eq (3.12) into Eq (3.10), we obtain

Using Cauchy-Schwarz inequality, we have

and

From Eq (3.14), Eq (3.15) and Eq (3.13), we obtain

Based on Lemma 3.3 and Lemma 3.4, we have

Applying Cauchy inequality, we get

Substituting Eq (3.18) into Eq (3.16) yields that

where Lemma 3.2 has been used. That is,

where $\mu = \Gamma (2 - {\alpha}){\tau ^{\alpha}}$. According to Lemma 3.1, we have

and

Substituting Eq (3.20) into Eq (3.19) gives that

Denote

Now, we prove by the mathematical induction method that

It holds obviously when $n = 0$. Assuming Eq (3.21) is valid for $n = 1, 2, \cdots, m-1$, then we have

This proof is completed.

Applying the Theorem 3.5, we can immediately obtain the stability result.

Theorem 3.6 (Stability). The difference scheme (2.10)–(2.13) is unconditionally stable with respect to the initial value $\varphi$ and the source term $f$.

Similarly, from Theorem 3.5, we can easily prove the solvability of the proposed scheme.

Theorem 3.7 (Solvability). The difference scheme (2.10)–(2.13) is uniquely solvable.

Proof. It suffices to prove the homogeneous linear system

has only a trivial solution. Applying Theorem 3.1, we have $||u^n||^2\le ||u^0||^2 = 0$. So $u_i^n\equiv0$ for $0\le i\le M$, which completes the proof.

Next, we focus on the convergence of the difference scheme (2.10)–(2.13). Denote

Theorem 3.8 (Convergence). Assume that $u(x, t)\in C^{6, 3}_{x, t}([0, L]\times[0, T])$ and $\{ u_i^n\}$ are solution of the problem (1.1)–(1.4) and the difference scheme Eqs (2.10)–(2.13) respectively. Then there exists a positive constant $C$ such that

Proof. From Eq (2.7), Eq (2.8) and Eqs (2.10)–(2.13), we have the error equations as

Applying Theorem 3.5, we get

Noticing Eq (2.9), we get

which shows that Eq (3.22) is valid with

This proof is completed.

4.   The $\mathcal{FL}2 - 1_{\sigma}$ scheme

Although the $\mathcal{L}2 - 1_{\sigma}$ scheme (2.10)–(2.13) has accuracy of second order in time, it is not conducive to calculation due to it needs all history data to get the solution at current time point. Also, here we present a fast scheme by applying the sum-of-exponentials approximation to the kernel function $t^{-\alpha}$.

The sum-of-exponentials approximation reads as:

Lemma 4.1 (^[22]). For the given $\alpha \in (0, 1)$, tolerance error $\varepsilon$, cut-off time step size $\tilde \tau$ and final time $T$, there are one positive integer $N_{exp}$, positive points $s_j$ and corresponding positive weights $w_j (j = 1, 2, \cdots, N_{exp})$ satisfying

and the number of exponentials needed is of the order

The fast evaluation of Caputo derivative, $\mathcal{FL}2 - 1_{\sigma}$ formula, is given as follows:

where $\lambda = \frac{{{\tau^{-\alpha}}}}{{{\Gamma(2-\alpha)}}}$, $\tilde w_j = \frac{{{1}}}{{{\Gamma(1-\alpha)}}}w_j$, and $\tilde V_j^{n}$ can be got form the following recursive relation

with $\tilde V_j^{0} = 0, (j = 1, 2, \cdots, N_{exp})$ and

The recursive relation (4.2) shows that the $\mathcal{FL}2 - 1_{\sigma}$ formula reduces significantly the computational complexity. Noticing that Eq (4.2) can be equivalently rewritten as the following summation form

thus we have

in which ${}^\mathcal{F} g_0^{(1, \alpha)} = \lambda a_0$, and for $n\ge 1$,

The equivalent expression (4.3) is more applicable in stability and convergence analysis.

With respect to the $\mathcal{FL}2 - 1_{\sigma}$ formula, we have the following some results.

Lemma 4.2 (^[22]). For any $\alpha \in (0, 1)$, and $u(t) \in C^{3}[0, T]$, it holds that

Lemma 4.3 (^[22]). Suppose ${\alpha} \in (0, 1), {}^\mathcal{F} g_k^{(n+1, \alpha)}$ is defined by Eq (4.4), then it holds that

Lemma 4.4 (^[22]). Suppose $u = \{ {u^n}\mid0 \le n \le N-1\}$ is a grid function defined on ${{\Omega}_\tau }$, then it holds that

Similar to the derivation of the $\mathcal{L}2 - 1_{\sigma}$ scheme (2.10)–(2.13), it follows from Eq (2.1), Eq (2.2) we have

and there exists a constant ${}^\mathcal{F} C_r$ such that

Omitting the small terms ${}^\mathcal{F} ({R_1})_i^n$ and ${}^\mathcal{F} ({R_2})_i^n$ in Eq (4.5) and Eq (4.6), from the boundary and initial conditions (2.3)–(2.4), we obtain the $\mathcal{FL}2 - 1_{\sigma}$ scheme for the problem (2.1)–(2.4) as follows

Next, we focus on the solvability, stability and convergence of the $\mathcal{FL}2 - 1_{\sigma}$ scheme.

Before the discussion, we first prove the following priori estimate.

Theorem 4.5. Suppose $\{w_i^n, z_i^n\mid0\le i\le M, 0\le n\le N\}$ satisfy the difference scheme

Then, we have

Proof. Similar to the proof of the Theorem 3.5, we can obtain from Eq (4.12) and Eq (4.13) that

Noticing that

we get

From Lemma 4.4, we can further obtain

Denote

Now, we prove by the mathematical induction that

It holds obviously when $n = 0$. Assuming Eq (4.19) is valid for $n = 1, 2, \cdots, m-1$, then we have

This proof is completed.

Based on Theorem 4.5, we can obtain the following stability theorems.

Theorem 4.6 (Stability). The $\mathcal{FL}2 - 1_{\sigma}$ scheme Eqs (4.8)–(4.11) is uniquely solvable, and unconditionally stable with respect to the initial value $\varphi$ and the source term $f$.

Theorem 4.7 (Convergence). Assume that $u(x, t)\in C^{6, 3}_{x, t}([0, L]\times[0, T])$ and $\{ u_i^n\}$ are solutions of the problem (1.1)–(1.4) and the $\mathcal{FL}2 - 1_{\sigma}$ scheme (4.8)–(4.11), respectively. Then there exists a positive constant $C$ such that

Proof. From Eq (2.7), Eq (2.8) and Eqs (4.8)–(4.11), we have the error equations as

Applying Theorem 4.5, we get

Noticing Eq (4.7), we get

which shows that Eq (4.20) is valid with $C = {{}^\mathcal{F} C}\sqrt {\frac{{{1}}}{{{}^\mathcal{F} C}}(\frac{{{L^4}}}{{18{C_1}}}+2{C_2})}$.

5.   Numerical experiments

5.1.   The nonuniform $L_1$ approximation

It should be pointed out that the proposed difference schemes are based on assumptions that the solution of problem is sufficiently smooth. But the singularity of the time fractional derivative may lead to weak singularity near the initial time which may influence the accuracy of numerical method. Thus, in order to overcome the possible singularity of the solution near $t = 0$, some related techniques have been developed, such as the initial correction techniques, non-uniform discretization and so on ^{[23,24,25,26]}. Because of this, a analogously scheme for the problem (1.1)–(1.4) based on the uniform mesh in space and graded mesh in time is first given as follows:

where

and

with $x_i = ih, t_n = (n/N)^rT, \tau_n = t_n-t_{n-1}$, where the constant mesh grading exponent $r \geq 1$. It should be noted that the graded mesh will be simplified to a uniform grid when $r = 1$.

5.2.   Numerical results

In this subsection, we rely on two numerical examples to verify the availability of the proposed methods.

Let

Example 5.1. First, we consider the following problem

where $\omega(x) = x^2+1$ and $f(x, t) = \cos(\pi x)\frac{{{\Gamma(4+\alpha)}}}{{{6}}}t^3+(t^{3+\alpha}+1)\big[\cos(\pi x)-2\pi ^2\cos(\pi x)+4x\pi^3\sin(\pi x)+(x^2+1)\pi^4\cos(\pi x)\big].$

It is not difficult to verify that the exact solutions of the problems 5.1 is $u(x, t) = \cos(\pi x)(t^{3+\alpha}+1)$, which satisfies the smoothness requirement in Theorems 3.8 and 4.7.

The numerical accuracy of both schemes are tested with respect to $\alpha = 0.25, 0.5, 0.75$, respectively. In calculation, we take $\varepsilon = 10^{-13}$, which is much less than $\tau^{2}$. The errors and convergence orders of the suggested two schemes are showed in Table 1. We can observe that the values of $\mbox{Ord}$ are always close to $2$, which means that the $\mathcal{L}2 - 1_{\sigma}$ scheme and the $\mathcal{FL}2 - 1_{\sigma}$ scheme have second order accuracy both in space and time for different $\alpha\in (0, 1)$. Table 2 lists the convergence orders of both schemes when $\tau = h$ and CPU time with $\alpha$ = 0.5. Obviously, the $\mathcal{FL}2 - 1_{\sigma}$ scheme is faster than the $\mathcal{L}2 - 1_{\sigma}$ scheme, especially for small $\tau$.

From the Tables 1, 2, we can see that these numerical results are consistent with the previous theoretical results. It shows the $\mathcal{L}2 - 1_{\sigma}$ scheme (2.10)–(2.13) and the $\mathcal{FL}2 - 1_{\sigma}$ scheme (4.8)–(4.11) are convergent with second order accuracy in space and time, and the $\mathcal{FL}2 - 1_{\sigma}$ scheme is more practical.

Example 5.2. Now, we consider the following problem

where $\kappa = 0, \omega(x) = e^{x}$ and

The exact solution of the example 5.2 is $u(x, t) = (t^{\alpha}+t^3)\sin x.$

The error and numerical accuracy of scheme (5.1)–(5.6) are listed in Tables 3–5 with respect to $\alpha = 0.4, 0.6, 0.8$ and some values of grading exponent $r$, respectively. We keep $M = 2N$ in calculation. These results show that the scheme (5.1)–(5.6) has accuracy of order $\alpha$ when $r = 1$, and accuracy of order $2-\alpha$ when $r \geq r_c = (2-\alpha)/\alpha$. The reason for this result is that the smoothness requirement of the solution in Theorems 3.8 and 4.7 is not satisfied.

Example 5.3. Finally, we consider the following space-time variable coefficient problem

where

The exact solution of above problem is also $u(x, t) = \cos(\pi x)(t^{3+\alpha}+1)$, while the variable coefficient function $\omega(x, t) = (xt)^2+1$ which depends on the variables $x$ and $t$.

Similar to the spatially variable coefficient problem, we apply the $\mathcal{L}2-{1_\sigma}$ scheme and the $\mathcal{FL}2-{1_\sigma}$ scheme to solve the problem in Example 5.3. Table 6 presents the numerical results. In calculation, we take $\varepsilon = 10^{-11}$. It is shown that the $\mathcal{L}2-{1_\sigma}$ scheme and the $\mathcal{FL}2-{1_\sigma}$ scheme are convergent with second order accuracy in space and time.

6.   Conclusion

In this paper, we propose two second order difference schemes in both space and time for solving the variable coefficient fourth-order fractional sub-diffusion equation subject to the second Dirichlet boundary conditions. The $\mathcal{L}2-{1_\sigma}$ formula and $\mathcal{FL}2-{1_\sigma}$ formula are applied to approximation the time Caputo fractional derivative. Compared with $\mathcal{L}2-{1_\sigma}$ scheme, the $\mathcal{FL}2-{1_\sigma}$ scheme shows the better computational efficiency. The unconditional stability, solvability and convergence of the two schemes are strictly proved by the discrete energy method. The nonuniform $L_1$ approximation for the such problem is also given. Numerical examples are given to verify the effectiveness of both schemes. It should be pointed out that the results in this paper can be directly extended to time-space variable coefficient problems if we constrain the coefficient function $\omega(w, t)$ satisfying that $0 < C_1\le \omega(w, t) \le C_2$.

Acknowledgments

This work described in this paper was supported by the Sichuan Science and Technology Program (Grant No. 2020YJ0110, Grant No. 2022JDTD0019), the National Natural Science Foundation of China (Grant No. 11801389) and the Laurent Mathematics Center of Sichuan Normal University and National-Local Joint Engineering Laboratory of System Credibility Automatic Verification (Grant No. ZD20220105).

Conflict of interest

The authors declare there is no conflict of interest.