齐次分数阶扩散方程的加权C-N格式及其修正

Weighted C-N Scheme of HomogeneousFractional Di?usion Equations and Its Correction

摘要: 对于齐次时间分数阶扩散方程，在精确解不光滑时，数值方法精度会下降。针对这种情况，本文提出加权Crank-Nicolson格式(简记为加权C-N格式)及其修正格式，在精确解不光滑时，修正原格式的第1步后，可恢复方法的时间2阶精度。本文接着给出详细的收敛性分析，并且数值算例验证了方法的有效性。

Abstract: For the homogeneous time fractional diffusion equation, the accuracy of the numerical method will decrease when the exact solution is not smooth enough. In this case, we consider a weighted Crank-Nicolson scheme (masked as weighted C-N scheme) and its correction. After correcting the first step of the weighted C-N scheme, the second-order time accuracy can be restored. Then we give the detailed convergence analysis, and numerical examples verify the effectiveness of the scheme.

1. 引言

分数阶微分方程的数值方法的精度常依赖于精确解的光滑性。针对精确解不光滑的情况，许多学者开始研究数值方法的修正格式，以保持离散格式的高精度。例如，Lubich [1] [2] [3] 给出基于1阶和2阶向后差分格式的两种修正方法，并给出收敛性分析。Yan [4] [5] 考虑L1格式的修正格式，并提出基于分段2次插值 的新型离散格式，然后给出其修正格式。Tadjeran [6] 提出分数阶扩散方程的2阶C-N格式，Jin [7] 接着考虑分数阶C-N格式的修正。在本文中，我们考虑齐次分数阶扩散方程的加权C-N格式，并针对精确解不光滑的情况提出比较简单的修正方法，只需修正原格式的第1步，即可保持格式的2阶时间精度。本文接着给出修正格式的收敛性分析，最后通过数值算例验证方法的有效性。

本文考虑如下齐次分数阶扩散方程：

$\{\begin{array}{l}{}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)-Au\left(x,t\right)=0,L<x<R,0<t\le T\\ u\left(x,0\right)=u_0,L<x<R\\ u\left(L,t\right)=b_L\left(t\right),u\left(R,t\right)=b_R\left(t\right),0<t\le T\end{array}$ (1)

其中 $0<\alpha \le 1$， ${}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)$ 是Caputo时间分数阶导数，其定义为：

${}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)=\{\begin{array}{l}\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-s\right)}^{-\alpha }{u}^{\prime }\left(s\right)\text{d}s,0<\alpha <1}\\ \frac{\text{d}u\left(x,t\right)}{\text{d}t},\alpha =1\end{array}$ (2)

算子A表示有界正则区域上的自伴正定二阶椭圆偏微分算子 [4]，满足

$\Vert {\left(z-A\right)}^{-1}\Vert \le C{\left|z\right|}^{-1},\forall z\in {\Sigma }_{\omega }=\left\{z\ne 0:\left|\arg z\right|<\omega \right\},\omega \in \left(\frac{\pi }{2},\pi \right)$ (3)

其中 $\Vert ·\Vert$ 表示 $L_2$ 范数，记 ${}_0{}^{R}D{}_t{}^{\alpha }u\left(x,t\right)$ 为Riemman-Liouville导数，则有 ${}_0{}^{R}D{}_t{}^{\alpha }\left(u\left(x,t\right)-u_0\right)={}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)$。

2. 加权C-N格式及其修正格式

令 $V\left(t\right)=u\left(t\right)-u_0$，则方程(1)可表示为：

${}_0{}^{C}D{}_t{}^{\alpha }V\left(x,t\right)-AV\left(x,t\right)=A{u}_0$ (4)

令 $\tau$ 为时间步长， $t_n=n\tau ,n=0,1,2,\cdots ,N,0\le t_n\le T$，h为空间步长， $x_i=L+ih,i=0,1,\cdots ,M,L\le x_i\le R$。 ${\bar{\partial }}_t^{\alpha }V\left(t_n\right):={\tau }^{-\alpha }{\displaystyle \underset{j=0}{\overset{n}{\sum }}{b}_{n-j}{V}^j}$ 表示 ${}_0{}^{R}D{}_t{}^{\alpha }u\left(x,t\right)$ 的向后Euler卷积逼近，其生成函数为：

$\tilde{b}\left(\xi \right)={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}{b}_j{\xi }^j}={\left(1-\xi \right)}^{\alpha }$ (5)

分数阶导数的逼近格式在 $t_n$ 时刻的时间精度为1阶，在 $t_n-\frac{\alpha }{2}\tau$ 时刻的时间精度为2阶，对 $AV\left(t\right)$ 在 $t_n$ 时刻和 $t_{n-1}$ 时刻作线性拉格朗日插值，可得 $t_n-\frac{\alpha }{2}\tau$ 时刻的逼近格式 $\left(1-\frac{\alpha }{2}\right)A{V}^n+\frac{\alpha }{2}A{V}^{n-1}$。带入方程(4)，可得加权C-N格式：

${\tau }^{-\alpha }{\displaystyle \underset{j=0}{\overset{n}{\sum }}{b}_{n-j}{V}^j}-\left(1-\frac{\alpha }{2}\right)A{V}^n-\frac{\alpha }{2}A{V}^{n-1}=A{u}_0,n=1,2,\cdots ,N$ (6)

在精确解不光滑的情况下，加权C-N格式达不到2阶时间精度。我们对加权C-N格式(6)的第1步得初值条件添加一个权系数 $c_0$，适当选取可使离散格式保持2阶精度，加权C-N修正格式如下：

$\{\begin{array}{l}{\tau }^{-\alpha }{\displaystyle \underset{j=0}{\overset{n}{\sum }}{b}_{n-j}{V}^j}-\left(1-\frac{\alpha }{2}\right)A{V}^n-\frac{\alpha }{2}A{V}^{n-1}=\left(1+c_0\right)A{u}_0,n=1\\ {\tau }^{-\alpha }{\displaystyle \underset{j=0}{\overset{n}{\sum }}{b}_{n-j}{V}^j}-\left(1-\frac{\alpha }{2}\right)A{V}^n-\frac{\alpha }{2}A{V}^{n-1}=A{u}_0,n=2,\cdots ,N\end{array}$ (7)

3. 加权C-N修正格式的收敛性分析

为了证明加权C-N修正格式的收敛性分析，我们先给出3个引理。

引理1：定义式子：

$\mu \left(z\right)=\frac{1-{\text{e}}^{-\tau z}}{{\left(1-\frac{\alpha }{2}+\frac{\alpha }{2}{\text{e}}^{-\tau z}\right)}^{\frac{1}{\alpha }}}\cdot \frac{c_0{\text{e}}^{-\tau z}+\frac{{\text{e}}^{-\tau z}}{1-{\text{e}}^{-\tau z}}}{1-\frac{\alpha }{2}+\frac{\alpha }{2}{\text{e}}^{-\tau z}}$ (8)

本文统一令 $z\in {\Gamma }_{\tau }$， ${\Gamma }_{\tau }=\left\{z\in \Gamma :\left|\Im z\right|\le \frac{\pi }{\tau }\right\}\subseteq {\Sigma }_{\omega }$ [1]，当权系数 $c_0=\frac{1-\alpha }{2}$ 时，则有：

$\left|\mu \left(z\right)-1\right|\le C{\tau }^2{\left|z\right|}^2$ (9)

其中C为正常数。

证明：令 $w=\tau z$， $f\left(w\right)=\mu \left(z\right)-1$，则 $\underset{x\to 0}{\lim }f\left(w\right)=0$， $\underset{x\to 0}{\lim }{f}^{\prime }\left(w\right)=c_0+\frac{\alpha -1}{2}=0$， $\underset{x\to 0}{\lim }{f}^{″}\left(w\right)\ne 0$，所以 $\mu \left(w\right)-1=O\left(w^2\right)$，即可证得 $\left|\mu \left(z\right)-1\right|\le C{\tau }^2{\left|z\right|}^2$。

引理2：定义式子：

$K\left(z\right)=z^{-1}{\left(z^{\alpha }-A\right)}^{-1}A$ (10)

$z_{\tau }=\frac{1-{\text{e}}^{-\tau z}}{{\left(1-\frac{\alpha }{2}+\frac{\alpha }{2}{\text{e}}^{-\tau z}\right)}^{\frac{1}{\alpha }}}$ (11)

则下式成立：

$\Vert K\left(z\right)\Vert \le C{\left|z\right|}^{-1}$ (12)

$\Vert K\left(z_{\tau }\right)\Vert \le C{\left|z\right|}^{-1}$ (13)

其中 $\Vert ·\Vert$ 表示 $L_2$ 范数，C为正常数。

证明：令 $w=\tau z$，则 $w\to 0$ 时， $\tau z_{\tau }=O\left(w\right)$，所以 $c\left|\tau z\right|\le \left|\tau z_{\tau }\right|\le C\left|\tau z\right|$，即 $c\left|z\right|\le \left|z_{\tau }\right|\le C\left|z\right|$。

由引理1可得 $z\in {\Sigma }_{\omega },z^{\alpha }\in {\Sigma }_{\omega }$，将 $z,z^{\alpha }$ 带入式(3)可得：

$\Vert {\left(z^{\alpha }-A\right)}^{-1}\Vert \le C{\left|z\right|}^{-\alpha }$ (14)

$\Vert {\left(z_{\tau }^{\alpha }-A\right)}^{-1}\Vert \le C{\left|z\right|}^{-\alpha }$ (15)

且下式成立：

${\left(z^{\alpha }-A\right)}^{-1}A=z^{\alpha }{\left(z^{\alpha }-A\right)}^{-1}-I$ (16)

带入式(10)即可证得式(12)和(13)：

$\Vert K\left(z\right)\Vert ={\left|z\right|}^{-1}\Vert z^{\alpha }{\left(z^{\alpha }-A\right)}^{-1}-I\Vert \le C{\left|z\right|}^{-1}$ (17)

$\Vert K\left(z_{\tau }\right)\Vert ={\left|z_{\tau }\right|}^{-1}\Vert z_{\tau }^{\alpha }{\left(z_{\tau }^{\alpha }-A\right)}^{-1}-I\Vert \le C{\left|z\right|}^{-1}$ (18)

引理3： $z,z^{\alpha }$ 分别由引理1和引理2定义，则下式成立：

$\Vert K\left(z_{\tau }\right)-K\left(z\right)\Vert \le C{\tau }^2\left|z\right|$ (19)

证明：令 $w=\tau z$，，则 $\underset{x\to 0}{\lim }g\left(w\right)=\underset{x\to 0}{\lim }{g}^{\prime }\left(w\right)=\underset{x\to 0}{\lim }{g}^{″}\left(w\right)=0$， $\underset{x\to 0}{\lim }{g}^{‴}\left(w\right)\ne 0$， $\left(0<\alpha <1\right)$，所以可得 $g\left(w\right)=O\left(w^3\right)$，即：

$\left|z_{\tau }-z\right|\le C{\tau }^2{z}^3$ (20)

由 $\left|z_{\tau }^{\alpha }-z^{\alpha }\right|\le C{\left|z\right|}^{\alpha -1}\left|z_{\tau }-z\right|\le C{\tau }^2{\left|z\right|}^{2+\alpha }$ ( [8] 引理B.1)，且由式(16)及可得：

$\begin{array}{l}\Vert K\left(z_{\tau }\right)-K\left(z\right)\Vert =\Vert z_{\tau }^{-1}{\left(z_{\tau }^{\alpha }-A\right)}^{-1}A-z^{-1}{\left(z^{\alpha }-A\right)}^{-1}A\Vert \\ \le {\left|z_{\tau }\right|}^{-1}\Vert {\left(z_{\tau }^{\alpha }-A\right)}^{-1}A-{\left(z^{\alpha }-A\right)}^{-1}A\Vert +\left|z_{\tau }^{-1}-z^{-1}\right|\Vert {\left(z^{\alpha }-A\right)}^{-1}A\Vert \\ \le {\left|z_{\tau }\right|}^{-1}\Vert z_{\tau }^{\alpha }{\left(z_{\tau }^{\alpha }-A\right)}^{-1}-z^{\alpha }{\left(z^{\alpha }-A\right)}^{-1}\Vert +\left|z_{\tau }^{-1}-z^{-1}\right|\Vert z^{\alpha }{\left(z^{\alpha }-A\right)}^{-1}-I\Vert \\ \le C{\left|z\right|}^{-1}\left|z_{\tau }^{\alpha }-z^{\alpha }\right|\Vert {\left(z_{\tau }^{\alpha }-A\right)}^{-1}\Vert +C{\left|z\right|}^{-1}{\left|z\right|}^{\alpha }\Vert {\left(z_{\tau }^{\alpha }-A\right)}^{-1}-{\left(z^{\alpha }-A\right)}^{-1}\Vert +C\frac{\left|z-z_{\tau }\right|}{\left|z{z}_{\tau }\right|}\\ \le C{\tau }^2{\left|z\right|}^{-1}{\left|z\right|}^{2+\alpha }{\left|z\right|}^{-\alpha }+C{\left|z\right|}^{\alpha -1}\Vert \left(z^{\alpha }-z_{\tau }^{\alpha }\right){\left(z_{\tau }^{\alpha }-A\right)}^{-1}{\left(z^{\alpha }-A\right)}^{-1}\Vert +C{\tau }^2\left|z\right|\le C{\tau }^2\left|z\right|\end{array}$ (21)

其中 ${\left(z_{\tau }^{\alpha }-A\right)}^{-1}-{\left(z^{\alpha }-A\right)}^{-1}=\left(z^{\alpha }-z_{\tau }^{\alpha }\right){\left(z_{\tau }^{\alpha }-A\right)}^{-1}{\left(z^{\alpha }-A\right)}^{-1}$，即证得引理3。

为了给出收敛性分析，接下来我们分别借助Laplace变换和Cauchy积分公式，给出齐次分数阶扩散方程(4)的精确解和加权C-N修正格式(7)的数值解。对方程(4)作Laplace变换可得：

$\hat{V}\left(z\right)=z^{-1}{\left(z^{\alpha }-A\right)}^{-1}A{u}_0$ (22)

对式(21)作Laplace逆变换可得方程(4)的精确解：

$V\left(t\right)=\frac{1}{2\pi i}{\displaystyle {\int }_{\Gamma }{\text{e}}^{tz}{z}^{-1}{\left(z^{\alpha }-A\right)}^{-1}A{u}_0\text{d}z}=\frac{1}{2\pi i}{\displaystyle {\int }_{\Gamma }{\text{e}}^{tz}K\left(z\right)u_0\text{d}z}$ (23)

其中 $\Gamma =\left\{z:\left|\arg z\right|=\theta \right\}\subseteq {\Sigma }_{\omega }$， $\theta \in \left(\frac{\alpha }{2},\frac{\alpha }{1+\alpha }\right)$ [4]。

考虑加权C-N修正格式(7)，式子两边同乘 ${\xi }^n$，并关于n求和， $n=1,2,\cdots$，令 $\tilde{V}\left(\xi \right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{V}^n}{\xi }^n$，于是 ${\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left({\displaystyle \underset{j=0}{\overset{n}{\sum }}{b}_{n-j}{V}^j}\right){\xi }^n}=\left({\displaystyle \underset{j=0}{\overset{\infty }{\sum }}{b}_j{\xi }^j}\right)\left({\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{V}^n{\xi }^n}\right)$，由式(5)可得：

${\tau }^{-\alpha }\tilde{b}\left(\xi \right)\tilde{V}\left(\xi \right)-\left(1-\frac{\alpha }{2}+\frac{\alpha }{2}\xi \right)A\tilde{V}\left(\xi \right)=\left(c_0\xi +\frac{\xi }{1-\xi }\right)A{u}_0$ (24)

$\tilde{V}\left(\xi \right)=\frac{c_0\xi +\frac{\xi }{1-\xi }}{1-\frac{\alpha }{2}+\frac{\alpha }{2}\xi }\cdot {\left[{\left(\frac{1-\xi }{\tau {\left(1-\frac{\alpha }{2}+\frac{\alpha }{2}\xi \right)}^{\frac{1}{\alpha }}}\right)}^{\alpha }-A\right]}^{-1}A{u}_0$ (25)

令 $\xi ={\text{e}}^{-\tau z}$，由Cauchy积分公式及式(8)、(10)和(11)可得加权C-N修正格式的数值解：

$\begin{array}{l}{V}^n=\frac{1}{2\pi i}{\displaystyle {\int }_{\left|\xi \right|=\rho }{\xi }^{-n-1}\tilde{V}\left(\xi \right)\text{d}\xi }=\frac{\tau }{2\pi i}{\displaystyle {\int }_{{\Gamma }_{\tau }}{\text{e}}^{t_{n}z}\tilde{V}\left({\text{e}}^{-\tau z}\right)\text{d}z}\\ =\frac{1}{2\pi i}{\displaystyle {\int }_{{\Gamma }_{\tau }}{\text{e}}^{t_{n}z}\mu \left(z\right)K\left(z_{\tau }\right)u_0\text{d}z}\end{array}$ (26)

定理1 $V\left(t_n\right),V^n$ 分别为 $t=t_n$ 时刻方程(4)的精确解和加权C-N修正格式(7)的数值解，则下式成立：

$\Vert V\left(t_n\right)-V^n\Vert \le C{\tau }^2{t}_n^{-2}\Vert u_0\Vert$ (27)

证明： $\Vert V\left(t_n\right)-V^n\Vert =\Vert \frac{1}{2\pi i}{\displaystyle {\int }_{{\Gamma }_{\tau }}{\text{e}}^{t_{n}z}\left[K\left(z\right)-\mu \left(z\right)K\left(z_{\tau }\right)\right]u_0\text{d}z}+\frac{1}{2\pi i}{\displaystyle {\int }_{\Gamma /{\Gamma }_{\tau }}{\text{e}}^{t_{n}z}K\left(z\right)u_0\text{d}z}\Vert \le \Vert I\Vert +\Vert II\Vert$，而 $\Vert \mu \left(z\right)K\left(z_{\tau }\right)-K\left(z\right)\Vert \le \left|\mu \left(z\right)-1\right|\Vert K\left(z_{\tau }\right)\Vert +\Vert K\left(z_{\tau }\right)-K\left(z\right)\Vert$，由引理1和引理2可得：

$\Vert \mu \left(z\right)K\left(z_{\tau }\right)-K\left(z\right)\Vert \le C{\tau }^2\left|z\right|$ (28)

接着考虑 $\Vert I\Vert ,\Vert II\Vert$，令 $r=t_n\left|z\right|$，c为正常数，则：

$\begin{array}{c}\Vert I\Vert \le \frac{1}{2\pi }{\displaystyle {\int }_{{\Gamma }_{\tau }}\left|{\text{e}}^{t_{n}z}\right|\Vert \mu \left(z\right)K\left(z_{\tau }\right)-K\left(z\right)\Vert \cdot \Vert u_0\Vert \text{d}z}\\ \le C{\displaystyle {\int }_{{\Gamma }_{\tau }}\left|{\text{e}}^{t_{n}z}\right|{\tau }^2\left|z\right|\Vert u_0\Vert \text{d}z}\\ \le C{\tau }^2{\displaystyle {\int }_0^{\infty }{\text{e}}^{-cr}}r{t}_n^{-1}{t}_n^{-1}\Vert u_0\Vert \text{d}r\\ \le C{\tau }^2{t}_n^{-2}\Vert u_0\Vert \end{array}$ (29)

$\begin{array}{c}\Vert I\Vert \le \frac{1}{2\pi }{\displaystyle {\int }_{\Gamma /{\Gamma }_{\tau }}\left|{\text{e}}^{t_{n}z}\right|\Vert K\left(z\right)\Vert \Vert u_0\Vert \text{d}z}\\ \le C{\displaystyle {\int }_{\frac{1}{\tau }}^{\infty }{\text{e}}^{-c{t}_n\left|z\right|}}{\left|z\right|}^{-2}\left|z\right|\text{d}z\\ \le C{\tau }^2{\displaystyle {\int }_0^{\infty }{\text{e}}^{-cr}}r{t}_n^{-1}{t}_n^{-1}\Vert u_0\Vert \text{d}r\\ \le C{\tau }^2{t}_n^{-2}\Vert u_0\Vert \end{array}$ (30)

则可证得 $$，即权系数取 $c_0=\frac{1-\alpha }{2}$ 时，加权C-N修正格式为时间2阶精度。

4. 数值算例

数值算例1：考虑如下齐次分数阶扩散方程：

$\{\begin{array}{l}{}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)-\frac{1}{4{\pi }^2}\frac{{\partial }^{2}u}{\partial t^2}=0,0<x<1,0<t\le T\\ u\left(x,0\right)=\sin \left(2\pi x\right),0<x<1\\ u\left(0,t\right)=u\left(1,t\right)=0,0<t\le T\end{array}$ (31)

该方程的精确解为 $u\left(x,t\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-t^{\alpha }\right)}^k}{\Gamma \left(\alpha k+1\right)}}\sin \left(2\pi x\right)$，u在 $t=0$ 时不光滑。不同时间步长的误差如表1所示，不同 $\alpha$ 下的误差与对应步长的对数关系如图1所示，其中空间步长取 $h={10}^{-4},T=1$，统一取离散误差为 $L_2$ 误差 $\Vert V\left(t_N\right)-V^N\Vert$。

Table 1. L 2 error of weighted C-N scheme (6) and weighted C-N modified scheme (7)

表1. 加权C-N格式(6)和加权C-N修正格式(7)的 $L_2$ 误差

表1给出不同 $\alpha$ 及不同时间步长取值下，两种方法所得的误差。方法(a)为加权C-N格式离散，方法(b)为加权C-N修正格式离散，由表1可以看出，两种格式误差均收敛，且修正后的误差更小，方法更精确。

Figure 1. The relationship between error and time step

图1. 误差与时间步长的关系

图1给出不同 $\alpha$ 下误差与时间步长的关系，两坐标均为对数坐标，其中图1(a)为加权C-N格式，图1(b)为加权C-N修正格式。正如理论证明的结果一样，加权C-N格式只有1阶时间精度，而修正后的格式可达到2阶时间精度。

数值算例2：考虑如下齐次分数阶扩散方程：

$\{\begin{array}{l}{}_0{}^{C}D{}_t{}^{\alpha }u\left(x,t\right)-\frac{1}{4{\pi }^2}\frac{{\partial }^{2}u}{\partial t^2}=0,0<x<1,0<t\le T\\ u\left(x,0\right)=x\left(1-x\right),0<x<1\\ u\left(0,t\right)=u\left(1,t\right)=0,0<t\le T\end{array}$ (32)

该方程的精确解为 $u\left(x,t\right)=\frac{8}{{\pi }^3}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{1}{{\left(2n+1\right)}^3}{E}_{\alpha }}\left(-{\left(2n+1\right)}^2{\pi }^2{t}^{\alpha }\right)\sin \left(\left(2n+1\right)\pi x\right)$，u在 $t=0$ 时不光滑。取 $\alpha =0.\text{1}$，不同时间步长下的两种格式的误差和精度如表2所示。

Table 2. Error and time accuracy of weighted C-N scheme (6) and weighted C-N modified scheme (7)

表2. 加权C-N格式(6)和加权C-N修正格式(7)的误差和时间精度

由表2可以看出，随着时间步长减小，两方法的误差均减小，但加权C-N格式只有1阶时间精度，而修正格式的精度可达到2阶。

基金项目

国家自然科学基金青年科学基金项目(61802129)，广东省国家青年基金纵向协同项目(2018A030310381)。