A fast method for variable-order space-fractional diffusion equations

TL;DR

This paper introduces a fast, efficient numerical method for solving variable-order space-fractional diffusion equations, overcoming structural challenges in the stiffness matrix and significantly reducing computational complexity.

Contribution

The paper presents a novel divided-and-conquer collocation method that approximates the stiffness matrix with a sum of Toeplitz matrices, improving computational efficiency for variable-order fractional PDEs.

Findings

Achieves $O(kN ext{log}^2 N)$ memory usage.

Achieves $O(k N ext{log}^3 N)$ computational complexity.

Demonstrates effectiveness through numerical experiments.

Abstract

We develop a fast divided-and-conquer indirect collocation method for the homogeneous Dirichlet boundary value problem of variable-order space-fractional diffusion equations. Due to the impact of the space-dependent variable order, the resulting stiffness matrix of the numerical approximation does not have a Toeplitz-like structure. In this paper we derive a fast approximation of the coefficient matrix by the means of a sum of Toeplitz matrices multiplied by diagonal matrices. We show that the approximation is asymptotically consistent with the original problem, which requires $O(kN\log^2 N)$ memory and $O(k N\log^3 N)$ computational complexity with $N$ and $k$ being the numbers of unknowns and the approximants, respectively. Numerical experiments are presented to demonstrate the effectiveness and the efficiency of the proposed method.

Tables

Equations

$$\begin{array}[]{c}-u^{\prime\prime}(x)-d\,{}_{0}^{C}D_{x}^{\alpha}u(x)=f(x),~{}~{}x\in[0,1];\\[3.61371pt] u(0)=u(1)=0,\end{array}$$

$${}_0^C{D}_x^{\alpha }g(x):=\frac{1}{\Gamma (2-\alpha )}{\int }_0^x\frac{g^{\prime \prime }(s)}{(x-s{)}^{\alpha -1}}ds.$$

$$\begin{array}[]{c}-u^{\prime\prime}(x)-d(x)~{}{}_{0}^{C}D_{x}^{\alpha(x)}u(x)=f(x),~{}~{}x\in(0,1);\\[7.22743pt] u(0)=u(1)=0,\end{array}$$

$${}_0{D}_x^{\alpha (x)}g(x):=\frac{1}{\Gamma (2-\alpha (x))}{\int }_0^x\frac{g^{\prime \prime }(s)}{(x-s{)}^{\alpha (x)-1}}ds$$

$$v(x)+\frac{d(x)}{\Gamma (2-\alpha (x))}{\int }_0^x\frac{v(s)}{(x-s{)}^{\alpha (x)-1}}ds=-f(x).$$

$$u(x)={\int }_0^{x}v(s)(x-s)ds-x{\int }_0^{1}v(s)(1-s)ds.$$

$$v_h(x)=\sum _{n=0}^{N+1}{v}_n{\varphi }_n(x),v_n:=v_h(x_n).$$

$$\begin{array}[]{l}\displaystyle v_{h}(x_{n})+\frac{d(x_{n})}{\Gamma(2-\alpha(x_{n}))}\int_{0}^{x_{n}}\frac{v_{h}(s)}{(x_{n}-s)^{\alpha(x_{n})-1}}ds=-f(x_{n});\end{array}$$

$$\begin{array}[]{l}\displaystyle u_{h}(x):=\int_{0}^{x}v_{h}(s)(x-s)ds-x\int_{0}^{1}v_{h}(s)(1-s)ds,\quad x\in[0,1].\end{array}$$

$$\begin{array}[]{l}\displaystyle\|u-u_{h}\|:=\max_{1\leq n\leq N}|u(x_{n})-u_{h}(x_{n})|\leq\left\{\begin{array}[]{ll}\displaystyle QN^{-(3-\alpha_{max})},&~{}~{}\alpha(0)>1,\\[5.42003pt] \displaystyle QN^{-2},&~{}~{}\alpha(0)=1.\end{array}\right.\end{array}$$

$$\begin{array}[]{l}\displaystyle\mathbf{A}\mathbf{v}=\mathbf{f},\end{array}$$

$$\begin{array}[]{l}\displaystyle A_{i,j}=\left\{\begin{array}[]{ll}\displaystyle 1+\frac{d(x_{i})h^{2-\alpha(x_{i})}}{\Gamma(4-\alpha(x_{i}))},&i=j,\\[10.84006pt] \displaystyle\frac{d(x_{i})h^{2-\alpha(x_{i})}}{\Gamma(4-\alpha(x_{i}))}t_{i,j},&1\leq j\leq i-1,\\[10.84006pt] \displaystyle 0,&\mbox{otherwise},\end{array}\right.\end{array}$$

$$t_{i,j}=(i-j-1{)}^{3-\alpha (x_i)}-2(i-j{)}^{3-\alpha (x_i)}+(i-j+1{)}^{3-\alpha (x_i)},$$

$$\begin{array}[]{l}\displaystyle f_{n}=-f(x_{n})-d(x_{n})v_{0}\Big{(}\frac{x_{n}^{2-\alpha(x_{n})}}{\Gamma(3-\alpha(x_{n}))}+\frac{(x_{n}-x_{1})^{3-\alpha(x_{n})}-x_{n}^{3-\alpha(x_{n})}}{h\Gamma(4-\alpha(x_{n}))}\Big{)}.\end{array}$$

$$\begin{array}[]{l}\displaystyle u_{1}=\frac{h^{2}(2v_{0}+v_{1})}{6}-hI,\\[7.22743pt] \displaystyle u_{n+1}=u_{n}+hV_{n}+\frac{h^{2}(2v_{n}+v_{n+1})}{6}-hI,~{}~{}1\leq n\leq N-1,\end{array}$$

$$\begin{array}[]{ll}\displaystyle V_{0}=\frac{h(v_{0}+v_{1})}{2},&\displaystyle V_{n}=V_{n-1}+\frac{h(v_{n}+v_{n+1})}{2},~{}~{}1\leq n\leq N-1.\end{array}$$

$$\begin{array}[]{rl}\displaystyle u_{n+1}&\displaystyle=\int_{0}^{x_{n+1}}v_{h}(s)(x_{n+1}-s)ds-x_{n+1}\int_{0}^{1}v_{h}(s)(1-s)ds\\[10.84006pt] &\displaystyle=\int_{0}^{x_{n}}v_{h}(s)(x_{n+1}-s)ds+\int_{x_{n}}^{x_{n+1}}v_{h}(s)(x_{n+1}-s)ds-x_{n+1}I\\[10.84006pt] &\displaystyle=\int_{0}^{x_{n}}v_{h}(s)(x_{n}-s)ds+h\int_{0}^{x_{n}}v_{h}(s)ds\\[10.84006pt] &\displaystyle\quad\quad+\frac{h^{2}(2v_{n}+v_{n+1})}{6}-x_{n+1}I.\end{array}$$

$$\begin{array}[]{rl}\displaystyle u_{n}&\displaystyle=\int_{0}^{x_{n}}v_{h}(s)(x_{n}-s)ds-x_{n}I.\end{array}$$

$$\begin{array}[]{ll}\displaystyle u_{n+1}&\displaystyle=u_{n}+x_{n}I+hV_{n}+\frac{h^{2}(2v_{n}+v_{n+1})}{6}-x_{n+1}I\\[9.03374pt] \displaystyle&\displaystyle=u_{n}+hV_{n}+\frac{h^{2}(2v_{n}+v_{n+1})}{6}-hI,\end{array}$$

$$\begin{array}[]{l}\displaystyle t_{i,j}\approx 2(i-j)^{3-\bar{\alpha}}\Big{[}1+(\bar{\alpha}-\alpha(x_{i}))\ln(i-j)+\frac{(\bar{\alpha}-\alpha(x_{i}))^{2}}{2!}\ln^{2}(i-j)\\[10.84006pt] \qquad\quad\displaystyle+\cdots+\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{s}}{s!}\ln^{s}(i-j)\Big{]}\cdot\bigg{[}{{3-\alpha(x_{i})}\choose{2}}\frac{1}{(i-j)^{2}}\\[10.84006pt] \displaystyle\qquad\quad+{{3-\alpha(x_{i})}\choose{4}}\frac{1}{(i-j)^{4}}+\cdots+{3-\alpha(x_{i})\choose 2k}\frac{1}{(i-j)^{2k}}\bigg{]}\end{array}$$

$$\begin{array}[]{l}\displaystyle R^{i,j}_{s,k}=(i-j)^{3-\bar{\alpha}}\hat{R}^{i,j}_{s+1}\bigg{[}\Big{(}1-\frac{1}{i-j}\Big{)}^{3-\alpha(x_{i})}-2+\Big{(}1+\frac{1}{i-j}\Big{)}^{3-\alpha(x_{i})}\,\bigg{]}\\[12.64746pt] \displaystyle\qquad\qquad\qquad+(i-j)^{3-\alpha(x_{i})}R^{i,j}_{2k+2}-(i-j)^{3-\bar{\alpha}}\hat{R}^{i,j}_{s+1}R^{i,j}_{2k+2}.\end{array}$$

$$\begin{array}[]{l}\displaystyle(i-j\pm 1)^{3-\alpha(x_{i})}=(i-j)^{3-\alpha(x_{i})}\Big{(}1\pm\frac{1}{i-j}\Big{)}^{3-\alpha(x_{i})}\\[10.84006pt] \displaystyle\quad=(i-j)^{3-\alpha(x_{i})}\bigg{[}1\pm{{3-\alpha(x_{i})}\choose{1}}\frac{1}{i-j}+{{3-\alpha(x_{i})}\choose{2}}\frac{1}{(i-j)^{2}}\\[10.84006pt] \displaystyle\qquad\pm{{3-\alpha(x_{i})}\choose{3}}\frac{1}{(i-j)^{3}}+\cdots+{3-\alpha(x_{i})\choose 2k}\frac{1}{{(i-j)}^{2k}}\\[10.84006pt] \displaystyle\qquad\pm{3-\alpha(x_{i})\choose 2k+1}\frac{1}{{(i-j)}^{2k+1}}+R_{2k+2}^{i,j,\pm}\bigg{]},\end{array}$$

$$R_{2k+2}^{i,j,\pm}:={3-\alpha(x_{i})\choose 2k+2}\Big{(}1\pm\frac{\theta_{\pm}}{i-j}\Big{)}^{1-\alpha(x_{i})-2k}\frac{1}{(i-j)^{2k+2}},~{}~{}\theta_{\pm}\in(0,1).$$

$$\begin{array}[]{l}\displaystyle t_{i,j}=(i-j)^{3-\alpha(x_{i})}\bigg{[}{{3-\alpha(x_{i})}\choose{2}}\frac{2}{(i-j)^{2}}+{{3-\alpha(x_{i})}\choose{4}}\frac{2}{(i-j)^{4}}+\cdots\\[10.84006pt] \displaystyle\qquad\qquad\qquad+{3-\alpha(x_{i})\choose 2k}\frac{2}{(i-j)^{2k}}+R^{i,j}_{2k+2}\bigg{]},\\[-18.06749pt] \end{array}$$

$$\begin{array}[]{ll}\displaystyle R^{i,j}_{2k+2}&:\displaystyle=R_{2k+2}^{i,j,-}+R_{2k+2}^{i,j,+}={3-\alpha(x_{i})\choose 2k+2}\frac{1}{(i-j)^{2k+2}}\\[14.45377pt] \displaystyle&\displaystyle\qquad\times\bigg{[}\Big{(}1-\frac{\theta_{-}}{i-j}\Big{)}^{1-\alpha(x_{i})-2k}+\Big{(}1+\frac{\theta_{+}}{i-j}\Big{)}^{1-\alpha(x_{i})-2k}\,\bigg{]}.\end{array}$$

$$\begin{array}[]{rl}\displaystyle(i-j)^{3-\alpha(x_{i})}&\displaystyle=(i-j)^{\bar{\alpha}-\alpha(x_{i})}(i-j)^{3-\bar{\alpha}}\\[1.8063pt] &\displaystyle=\bigg{[}1+\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}\ln(i-j)+\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{2}}{2!}\ln^{2}(i-j)\\[7.22743pt] &\displaystyle\quad+\cdots+\frac{(\bar{\alpha}-\alpha(x_{i}))^{s}}{s!}\ln^{s}(i-j)+\hat{R}^{i,j}_{s+1}\bigg{]}(i-j)^{3-\bar{\alpha}},\end{array}$$

$$\begin{array}[]{l}\displaystyle\hat{R}^{i,j}_{s+1}=\frac{(\bar{\alpha}-\eta)^{s+1}\ln^{s+1}(i-j)}{(s+1)!},\quad\eta\in\left\{\begin{array}[]{l}\displaystyle(\bar{\alpha},\alpha(x_{i})),~{}~{}\mbox{if}~{}\bar{\alpha}\leq\alpha(x_{i}),\\[7.22743pt] (\alpha(x_{i}),\bar{\alpha}),~{}~{}\mbox{otherwise}.\end{array}\right.\end{array}$$

$$\begin{array}[]{l}\displaystyle t_{i,j}=2(i-j)^{3-\bar{\alpha}}\bigg{[}1+\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}\ln(i-j)+\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{2}}{2!}\ln^{2}(i-j)\\[10.84006pt] \qquad\quad\displaystyle+\cdots+\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{s}}{s!}\ln^{s}(i-j)+\hat{R}_{s+1}^{i,j}\bigg{]}\bigg{[}{{3-\alpha(x_{i})}\choose{2}}\frac{1}{(i-j)^{2}}\\[10.84006pt] \displaystyle\qquad\quad+{{3-\alpha(x_{i})}\choose{4}}\frac{1}{(i-j)^{4}}+\cdots+{3-\alpha(x_{i})\choose 2k}\frac{1}{(i-j)^{2k}}+R^{i,j}_{2k+2}\bigg{]}\\[10.84006pt] \displaystyle\quad\,\,\,=2(i-j)^{3-\bar{\alpha}}\bigg{[}1+(\bar{\alpha}-\alpha(x_{i}))\ln(i-j)+\frac{(\bar{\alpha}-\alpha(x_{i}))^{2}}{2!}\ln^{2}(i-j)\\[10.84006pt] \qquad\quad\displaystyle+\cdots+\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{s}}{s!}\ln^{s}(i-j)\bigg{]}\bigg{[}{{3-\alpha(x_{i})}\choose{2}}\frac{1}{(i-j)^{2}}\\[10.84006pt] \displaystyle\qquad\quad+{{3-\alpha(x_{i})}\choose{4}}\frac{1}{(i-j)^{4}}+\cdots+{3-\alpha(x_{i})\choose 2k}\frac{1}{(i-j)^{2k}}\bigg{]}+R^{i,j}_{s,k},\end{array}$$

$$\begin{array}[]{l}\displaystyle R_{s,k}^{i,j}=2(i-j)^{3-\bar{\alpha}}\bigg{[}\hat{R}^{i,j}_{s+1}\bigg{(}R^{i,j}_{2k+2}+\sum_{m=0}^{k}{{3-\alpha(x_{i})}\choose{2m}}\frac{1}{(i-j)^{2m}}\bigg{)}\\[10.84006pt] \displaystyle\qquad\qquad+R^{i,j}_{2k+2}\bigg{(}\sum_{m=0}^{s}\frac{\big{(}\bar{\alpha}-\alpha(x_{i})\big{)}^{m}}{m!}\ln^{m}(i-j)+\hat{R}_{s+1}\bigg{)}-R^{i,j}_{2k+2}\hat{R}^{i,j}_{s+1}\bigg{]}.\end{array}$$

$$\tilde{\mathbf{A}}\mathbf{v}=\mathbf{f},\tilde{\mathbf{A}}:=({\tilde{A}}_{i,j}{)}_{i,j=1}^{N+1}.$$

Full text

1 Introduction

Field tests showed that space-fractional diffusion equations (sFDEs) provide more accurate descriptions of challenging phenomena of superdiffusive transport and long range interaction, which occur in solute transport in heterogeneous porous media and other applications, than integer-order diffusion equations do BenSchMee ; DelCar ; MetKla04 ; SchBenMee01 . In fact, integer-order diffusion equations were derived if the underlying independent and identically distributed particle movements have (i) a mean free path and (ii) a mean waiting time. In this case, the central limit theorem concludes that the (normalized) partial sum of the independent and identically distributed particle movements converges to Brownian motions. The probability density distribution of finding a particle somewhere in space is Gaussian, which satisfies the classical Fickian diffusion equation MeeSik ; MetKla04 .

Note that assumptions (i) and (ii) hold for diffusive transport of solute in homogeneous porous media, where solute plumes were observed to decay exponentially Bear61 ; Bear72 and so can be described accurately by integer-order diffusion equations. However, field tests showed that solute transport in heterogeneous aquifers often exhibit highly skewed and power-law decaying behavior, while sFDEs were derived under the assumption that the solutions have such behavior BenSchMee ; MeeSik ; MetKla04 . This is why sFDEs can accurately describe the solute transport in heterogeneous media more accurately than integer-order diffusion equations do. Consequently, they have attracted extensive research activities in the last few decades ErvRoo05 ; LiZha ; LiChe ; LiuAnh .

However, sFDEs present new mathematical and numerical issues that are not common in the context of integer-order diffusion equations. Because of their nonlocal nature, numerical discretizations of sFDEs usually yield dense or full stiffness matrices Den ; LiZha ; LiuAnh ; Roo . A direct solver typically has $O(N^{3})$ computational complexity and $O(N^{2})$ memory requirement. A conventional Krylov subspace iterative method has $O(N^{2})$ computational complexity per iteration, but may diverge due to significant amount of round-off errors WanDu13a . In any case, the significantly increased computational complexity of numerical discretizations of sFDEs compared to their integer-order analogues is deemed computationally intractable for realistic simulations in multiple space dimensions, especially when parameter learning or control of the systems is involved.

It was discovered that the stiffness matrices of the numerical discretizations of constant-order sFDEs on a uniform partition typically possess a Toeplitz-like structure WanWanSir , which reduces the memory requirement from $O(N^{2})$ to $O(N)$ and computational complexity from $O(N^{3})$ to $O(N\log N)$ per Krylov subspace iteration via the discrete fast Fourier transform. Furthermore, different preconditioners were employed to futher improve the computational efficiency and even convergence behavior BaiPan ; Daniele ; JinF ; Jin3 ; LinNGSun ; LinNGSun1 ; PanNgWang ; WanDu13a ; ZJinL .

However, Ervin et al. ErvHeu16 proved that the one-dimensional constant-order constant-coefficient linear sFDEs with smooth right-hand side generate solutions with singularity at the end points of the spatial interval, which is in sharp contrast to their integer-order analogues and makes the error estimates of their numerical approximations derived under full regularity assumptions inappropriate. The singularity of the solutions to constant-order sFDEs seems to be physically irrelevant to the diffusive process of the solute, and occurs due to the incompatibility between the nonlocality of the power law decaying tails of the sFDEs inside the domain and the locality of the imposed classical boundary conditions. Intuitively, a physically relevant sFDE model should not only properly model the anomalous transport of solutes in heterogeneous porous media, but also correct the non-physical behavior of solutions to the existing sFDEs and thus maintain the smoothing nature of the diffusive transport process.

Recently, the wellposedness and smoothing properties of a variable-order linear sFDE was analyzed in ZheWan : If the variable order has an integer limit at the boundary then the solutions have the full regularity as their integer-order analogues do; Otherwise, the solutions exhibit certain singularity at the boundary as their constant-order sFDE analogues do. Thus, the variable-order sFDE provides a feasible approach to resolve the non-physical singularity of solutions to constant-order sFDEs near the boundary while retaining their advantages. In fact, variable-order sFDEs have been used in many applications SunChaZha ; SunCheChe , as the variable order is closely related to the fractal dimension of the porous media via the Hurst index EmbMae ; MeeSik and so can account for the changes of the geometrical structure or properties of the media.

Some numerical studies of variable-order FDEs can be found in the literature in recent years. First-order convergence rates were proved for finite difference methods for space-time-dependent variable-order space-fractional advection-diffusion equations in one space dimension under (the artificially assumed) full regularity assumptions of the true solutions without addressing the singularity issue of the problem ZhuLiu . A spectral collocation method using weighted Jacobi polynomials was derived for variable-order sFDEs ZenZhaKar , in which numerical experiments were presented to demonstrate the utility of the proposed method. The stability and convergence of an implicit alternating direct method was proved in CheLiuBur under the assumptions of the smoothness of the true solutions and (somewhat artificial) monotonicity of the diffusivity coefficients. However, due to the impact of variable order of the FDEs, the numerical discretizations of variable-order sFDEs no longer have Toeplitz-like stiffness matrices, so the fast solvers developed for constant-order sFDEs do not apply. Also, the spectral methods do not have diagonal stiffness matrices.

In this paper we develop a fast numerical solution technique for an indirect collocation method to the homogeneous Dirichlet boundary-value problem of a one-sided variable-order sFDE in one space dimension. We approximate the stiffness matrix by a finite sum of Toeplitz-like matrices that is asymptotically convergent to the stiffness matrix. Then we develop a fast divided and conquer (DAC) solver for the approximated system by employing the Toeplitz-like structures of each summand to reduce the computational complexity from $O(N^{2})$ to $O(kN\log^{3}N)$ and the memory requirement from $O(N^{2})$ to $O(kN\log^{2}N)$.

The rest of the paper is organized as follows. In Section 2 we present the model problem and its numerical discretization. In Section 3 we approximate the coefficient matrix by a sum of Toeplitz-like matrices and analyze its asymptotic consistency. In Section 4 we develop a fast DAC method for the approximated system. We perform numerical experiments to test the performance of the method in the last section.

2 A variable-order sFDE model and its indirect collocation method

2.1 Model problem

An sFDE of order $1<\alpha<2$ was proposed in del to model the anomalously superdiffusive transport of solute in heterogeneous porous media

[TABLE]

where $u^{\prime\prime}$ refers to the second-order derivative of $u$ and the Caputo fractional derivative ${}_{0}^{C}D_{x}^{\alpha}g$ is defined by Pod

[TABLE]

Equation (1) is proposed based on the fact that a large amount of solute particles may travel through high permeability zones in a superdiffusive manner BenSchMee ; MeeSik , which may deviate from the transport of the solute particles in the bulk fluid phase that undergo a Fickian diffusive transport Bear72 . Therefore, in model (1), the $-u^{\prime\prime}$ term represents the $1/(1+d)$ portion of the total solute mass undergoing the Fickian diffusive while the $-d\;{}_{0}^{C}D_{x}^{\alpha}u$ term refers to the $d/(1+d)$ portion of the total solute mass undergoing the superdiffusive transport in high permeability zones.

Note that in realistic applications, the reservoir may consist of different types of porous media that have different fractional dimensions. Hence, in this paper we consider the homogeneous Dirichlet boundary-value problem of the following one-sided variable-order linear sFDE as a variable extension of (1)

[TABLE]

where $1\leq\alpha_{min}\leq\alpha(x)\leq\alpha_{max}<2$ ($\alpha(x)\not\equiv\alpha_{min}$ if $\alpha_{min}=1$) and $d(x)\geq 0$ is the fractional diffusivity. The variable-order Caputo fractional derivative ${}_{0}D_{x}^{\alpha(x)}g$ is defined by ZenZhaKar ; ZhuLiu

[TABLE]

where $\Gamma(\cdot)$ refers to the Gamma function.

2.2 An indirect collocation method

We rewrite (2) in terms of $v(x):=u^{\prime\prime}(x)$

[TABLE]

Then the solution $u$ to model (2) can be obtained by postprocessing

[TABLE]

Let $0=x_{0}<x_{1}<\cdots<x_{N+1}=1$ be a uniform partition of $[0,1]$ with $x_{n}=nh$ for $n=0,1,\cdots,N+1$ and $h=1/(N+1)$. Let $\phi_{n}(x)$ be the piecewise-linear basis functions with $\phi_{n}(x_{n})=1$ and $\phi_{n}(x_{m})=0$ for $m\neq n$. Each element $v_{h}(x)$ in the space $S_{h}$ of continuous and piecewise linear functions on $x\in[0,1]$ can be represented by

[TABLE]

Then an indirect collocation method for model (2) reads:

Step 1

Find $v_{h}(x)\in S_{h}$ such that for $0\leq n\leq N+1$

[TABLE]

Step 2

Obtain an approximation $u_{h}(x)$ of $u(x)$ by

[TABLE]

The following error estimate of the method was proved ZheWan under the assumptions that $d$, $f$ and $\alpha$ are second-order continuously differentiable on $[0,1]$, and $1\leq\alpha_{min}\leq\alpha(x)\leq\alpha_{max}<2$ with $\alpha(x)\not\equiv\alpha_{min}$ if $\alpha_{min}=1$

[TABLE]

Instead, for $u$ sufficiently smooth, we have $\|u-u_{h}\|\leq QN^{-2}$.

2.3 Solution of the numerical scheme

For $n=0$ we obtain from (6) that $v_{0}=-f(0)$. We then plug (5) into (6) for $1\leq n\leq N+1$ and move the terms containing $v_{0}$ to the right-hand side of the resulting equation to get a linear system with the unknowns $\{v_{n}\}_{n=1}^{N+1}$

[TABLE]

where $\mathbf{A}=(A_{i,j})_{i,j=1}^{N+1}\in\mathbb{R}^{(N+1)\times(N+1)}$ is a lower triangular matrix of the form

[TABLE]

[TABLE]

$\mathbf{v}:=[v_{1},v_{2},\cdots,v_{N+1}]^{T}$ and $\mathbf{f}:=[f_{1},f_{2},\cdots,f_{N+1}]^{T}$ with $f_{n}$ given by

[TABLE]

With $v_{h}(x)$ obtained from (9), we plug $v_{h}(x)$ into (7) to obtain the approximation $u_{h}(x)$ to $u(x)$. In particular, the discrete approximation $\mathbf{u}=[u_{1},u_{2},\cdots,u_{N}]^{T}$ with $u_{n}:=u_{h}(x_{n})$ can be evaluated inductively as follows.

Theorem 2.1

The $\{u_{n}\}_{n=1}^{N}$ can be computed by

[TABLE]

where $I=\int_{0}^{1}v_{h}(s)(1-s)ds$ and $\{V_{n}\}_{n=0}^{N-1}$ is inductively generated by

[TABLE]

Proof

For (7) with $x=x_{1}$, we obtain the first equation of (12) by a direct calculation. For (7) with $x=x_{n+1}$ and $1\leq n\leq N-1$, we split the first integral on $(0,x_{n+1})$ to those on $(0,x_{n})$ and $(x_{n},x_{n+1})$ and split the kernel $(x_{n+1}-s)$ to $(x_{n}-s)+h$ in the first resulting integral to obtain

[TABLE]

Note that as $v_{h}(x)$ is a piecewise linear function, $V_{n}$ defined by (13) represents the integration of $v_{h}(x)$ over $(0,x_{n})$. Therefore, the second term on the right-hand side of (14) equals to $hV_{n}$. From (7) with $x=x_{n}$ we have

[TABLE]

We plug this into (14) to get

[TABLE]

which finishes the proof.

3 An approximated collocation scheme

As $\mathbf{A}$ is a lower triangular matrix, a direct solution of the numerical scheme requires $O(N^{2})$ storage and has $O(N^{2})$ computational complexity. Once we obtained $\mathbf{v}$, only $O(N)$ storage and $O(N)$ operations are needed to compute $\mathbf{u}$ by Theorem 2.1. That is, the main task lies in reducing the memory requirement and the computational complexity of (9). To do so, we approximate $\mathbf{A}$ by a finite sum of Toeplitz-like matrices.

Theorem 3.1

For $1\leq j\leq i-1$, $t_{i,j}$ defined by (11) can be approximated by

[TABLE]

for some $k,s\in\mathbb{N}^{+}$ with the residue $R^{i,j}_{s,k}$ given by

[TABLE]

Proof

We decouple the nonlinear dependence of $i$ and $j$ in $t_{i,j}$ (cf (11)) for $j\leq i-1$ by applying the Taylor’s expansion of power functions to the first and the third terms on the right-hand side of (11) as follows

[TABLE]

where $R_{2k+2}^{i,j,pm}$ are the remainders of the Taylor’s expansion

[TABLE]

Here $\theta_{\pm}=\theta_{\pm}(i,j,k)$. We plug (17) into (11) to obtain

[TABLE]

[TABLE]

Due to the impact of the variable order $\alpha(x)$, the matrix of entries $(i-j)^{3-\alpha(x_{i})}$ is not Toeplitz. The assumption on $\alpha(x)$ (preceding (8)) implies $\bar{\alpha}:=(\alpha_{max}+\alpha_{min})/2>1$. We use the Taylor’s expansion of $a^{x}$ with $a>0$ and $0<x<1$ to obtain

[TABLE]

[TABLE]

Substituting (21) into (19) yields

[TABLE]

[TABLE]

Then we replace the terms in the first and second brackets by (19) and (21), respectively, to obtain (16). Dropping the truncation error $R^{i,j}_{s,k}$ in (23), we get the approximation (15) of $t_{i,j}$ for $j\leq i-1$.

Replacing $t_{i,j}$ in the entries of $\mathbf{A}$ by the right-hand side of (15) leads to the corresponding approximated system

[TABLE]

Theorem 3.2

For $j\leq i-1$ and $k,s\in\mathbb{N}^{+}$ with $s>e\ln N/2-1$, the local truncation error $R^{i,j}_{s,k}$ can be bounded by

[TABLE]

Consequently, we have for $i-j\geq 2$

[TABLE]

That is, the approximation scheme (24) is asymptotically consistent with the original problem (9) with respect to $s$ and $k$.

Proof

We bound $R^{i,j}_{2k+2}$ and $\hat{R}^{i,j}_{s+1}$ respectively. The binomial coefficients in (20) can be represented in terms of the gamma functions as follows

[TABLE]

We use the asymptotic expansions of Gamma functions ((KilSri, , Eq. 1.5.15))

[TABLE]

to get

[TABLE]

We note that for $i-j=1$ the left-hand side of (17) with the minus sign in “$\pm$” vanishes, so $R^{i,j,-}_{2k+2}$, and thus the first term on the right-hand side of the equal sign“=” in (28), vanishes. We thus bound the remaining factors on the right hand of (20) by

[TABLE]

We incorporate the proceeding estimates to obtain the estimates of the local truncation error $R_{2k+2}^{i,j}$

[TABLE]

We note that $\eta$ in (22) satisfies $(\bar{\alpha}-\eta)^{s+1}\leq 2^{-(s+1)}$ and use the Stirling’s formula for gamma functions to get

[TABLE]

We thus bound $\hat{R}_{s+1}^{i,j}$ by

[TABLE]

We combine (29) and (30) to bound $R^{i,j}_{s,k}$ in (16) by

[TABLE]

where in the first inequality we have used the Taylor’s expansion as in (17) and the expression (18) for $R^{i,j,\pm}_{2k+2}$ with $k=0$

[TABLE]

By setting $\displaystyle s>\frac{e\ln N}{2}-1$ in the first term on the right-hand side of (31) we obtain (25) and (26).

4 A fast divided-and-conquer solver

We develop a fast DAC solver for the collocation system (24). Due to the impact of the variable order in the sFDE (2), the stiffness matrix $\mathbf{A}$ in (9) no longer has Toeplitz structure so the DAC algorithm MSun ; FMWang developed for constant-order FDEs cannot apply. We develop a fast DAC method by expressing the stiffness matrix $\mathbf{\tilde{A}}$ in (24) as

[TABLE]

By Theorem 3.2 the truncation errors $R^{i,j}_{s,k}$ are not necessarily small when $i-j$ is small, so the approximations of the corresponding entries may lose accuracy. To ensure the accuracy of the approximations, we evaluate entries on the right-up bands of $\mathbf{\Gamma}_{N^{\prime}}$ with width ${c=\lceil\log N\rceil}$ (i.e., colored entries in the matrix below) exactly as follows

[TABLE]

Thus, the matrix $\mathbf{\Psi}_{N^{\prime}}$ has $\frac{c(c+1)}{2}$ ($c=\lfloor\log N\rfloor$) nonzero entries totally, which indicates the matrix-vector multiplication $\mathbf{\Psi}_{N^{\prime}}\mathbf{v}_{N^{\prime}}$ can be computed exactly in $\frac{1}{2}\log^{2}N$ of memory and $\frac{1}{2}\log^{2}N$ of computational work.

Theorem 4.1

The sub-matrix $\tilde{\mathbf{\Gamma}}_{N^{\prime}}$ in (32) can be expressed as a sum of Toeplitz matrices multipied by diagonal matrices, so the corresponding system can be stored in $O(kN\log^{2}N)$ and can be solved in $O(kN\log^{3}N)$ operations by the fast approximated DAC method (see Algorithm 2).

Proof

We combine (33) with (15) to express $\mathbf{\hat{\Gamma}}_{N^{\prime}}$ as

[TABLE]

where $\mathbf{K}^{p,q}:=(K^{p,q}_{i})_{i=1}^{N^{\prime}}$ and $\mathbf{T}^{p,q}$ for $0\leq p\leq s,1\leq q\leq k$ are given by

[TABLE]

and

[TABLE]

with $\mathbf{t}_{c}^{p,q}:=(t_{ci}^{p,q})_{i=1}^{N^{\prime}}$ and $\mathbf{t}_{r}^{p,q}:=(t_{ri}^{p,q})_{i=1}^{N^{\prime}}$ the first column and the first row of $\mathbf{T}^{p,q}$, respectively

[TABLE]

It is well known that $\mathbf{T}^{p,q}\mathbf{v}_{N^{\prime}}$ can be performed in $O(N^{\prime}\log N^{\prime})$ operations via the discrete fast Fourier transform (FFT), which implies that $\mathbf{\hat{\Gamma}}_{N^{\prime}}\mathbf{v}_{N^{\prime}}$ can be evaluated in $O(ksN^{\prime}\log N^{\prime})$ operations. To store $\mathbf{\hat{\Gamma}}_{N^{\prime}}$, we only need to store $\mathbf{t}_{c}^{p,q}$, $\mathbf{t}_{r}^{p,q}$, and $\mathbf{K}^{p,q}$ that require $O(ksN^{\prime})$ memory. If we take $s=\lceil{\frac{e}{2}\ln N-1}\rceil=O(\log N)$ according to Theorem 3.2, the computational cost and the memory requirement of $\mathbf{\hat{\Gamma}}_{N^{\prime}}\mathbf{v}_{N^{\prime}}$ become $O(kN^{\prime}\log N\log N^{\prime})$ and $O(kN^{\prime}\log N)$, respectively. Then the total number of computations ${\tilde{\Theta}}_{N}$ and the storage $\tilde{\Phi}_{N}$ of the proposed fast method can be evaluated by

[TABLE]

By (33), the storage and the computational cost of $\mathbf{\Psi}_{N^{\prime}}\mathbf{v}_{N^{\prime}}$ are both $\frac{1}{2}\log^{2}N$. Thus the total computational cost $\epsilon_{\theta}$ and the storage $\epsilon_{\phi}$ are

[TABLE]

Then totally flops $\Theta_{N}$ and storage $\Phi_{N}$ of the fast approximated algorithm can be estimated by

[TABLE]

which finishes the proof.

Remark 1

It takes $O(N^{2})$ computational works to generate the entries of the original coefficient matrix $\mathbf{A}$, while in the fast approximated DAC method only $O(ks\log N)$ operations are needed to generate all the entries of the approximated coefficient matrix $\tilde{\mathbf{A}}$.

5 Numerical experiments

We perform numerical experiments to investigate the performance of the fast DAC method (FDAC) by comparing it with the Forward Substitution method (FS). The approximated system (24) will be solved and the convergence rates of the indirect collocation method, the CPU times (in seconds) of generating the coefficient matrix ($CPU_{M}$) and of solving the lower triangular linear system ($CPU_{S}$) will be recorded. We set the parameter $k$ appeared in Theorem 3.1 equal to $2$ throughout the experiments.

Experiment 1.

We set $d(x)=1$ and

[TABLE]

with $\alpha_{0}=1.2$, $\alpha_{1}=1.6$. The exact solution $u(x)=x^{4}(1-x)$ and the corresponding right hand term $f(x)$ is given by

[TABLE]

The results are presented in Table 1 and Table 2, from which we notice that

• •

The CPU time consumed by FS increases at about a quadruple rate between two consecutive numbers of collocation grids while the increment of the CPU time of the FDAC is almost liner;

• •

It is more efficient in the FDAC to generate the entries of the coefficient matrix than that in the FS. For instance, when $N=2^{15}$, the FS takes 352 seconds to compute the entries of $\mathbf{A}$ while the FDAC only requires $3.78$ seconds;

• •

when $N\geq 2^{11}$, the FDAC is more efficient than FS for solving the linear systems. For instance, for the case of $N=2^{15}$, the $CPU_{S}$ of FS is $99$ seconds while in the FDAC it is $1.5$ seconds;

• •

The FS is out of memory for $N\geq 2^{16}$ while the FDAC still works even for $N\geq 2^{19}$;

• •

The FDAC has almost the same accuracy and convergence rates as FS for relatively small $N$ while for large $N$ the convergence rate is affected by the round-off errors.

From the observations mentioned above, the presented FDAC has shown strong potentials for efficiently and effectively solving the variable-order space-fractional diffusion equations by reducing the memory requirement and improving the efficiency of generating entries of the coefficient matrix and solving the linear systems. This implies that the proposed fast method is particularly suitable for the large-scale simulations.

Experiment 2.

Let $d(x)=1$ and $\alpha(x)=(\alpha_{1}-\alpha_{0})x+\alpha_{0}$. We choose a solution with a boundary layer at $x=0$

[TABLE]

The corresponding right-hand side is

[TABLE]

Numerical results are presented in Table 3–5, which is strongly consistent with the observations in Experiment 1.

Acknowledgments

This work was funded by the OSD/ARO MURI Grant W911NF-15-1-0562, by the National Science Foundation under grants DMS-1620194 and by Natural Science Foundation of Shandong Province under grants ZR2019BA026.