[A1] 对一类分数阶扩散方程的离散化

\quad\quad现在对如下的一类扩散方程进行离散化:\begin{array}{l}u(x, y, t)_{t}=-\mathrm{div} \left(\mathbf{q}(x, y, t)\right) \text { in } \Omega \times(0, T), \\u(x, y, t)=0 \text { on } \partial \Omega \times(0, T), \\u(x, y, 0)=f(x, y) \text { in } \Omega,\end{array} 其中 \Omega=[a,b]\times[c,d] \subset \mathbb{R}^2,扩散通量 \mathbf{q}(x, y, t)=(q_x,q_y) 定义如下:\begin{array}{l}q_{x}=-c\left(u,\left|D_{x+}^{\alpha}u\right|\right)D_{x+}^{\alpha}u+c\left(u,\left|D_{x-}^{\alpha} u\right|\right) D_{x-}^{\alpha} u, \\[1ex]q_{y}=-c\left(u,\left|D_{y+}^{\alpha} u\right|\right)D_{y+}^{\alpha}u+c\left(u,\left|D_{y-}^{\alpha}u\right|\right)D_{y-}^{\alpha} u,\end{array} 其中 0<a\leqslant3D_{x+}^{\alpha}, D_{x-}^{\alpha},D_{y+}^{\alpha},D_{y-}^{\alpha} 是两个方向上的左、右 Riemann-Liouville 导数,定义如下:\begin{aligned} D_{x+}^{\alpha} u(x, y, t) &= \displaystyle \frac{1}{\varGamma(m-\alpha)}\left(\frac{\partial}{\partial x}\right)^{m} \int_{a}^{x} \frac{u(\xi, y, t)}{(x-\xi)^{\alpha+1-m}} d \xi \\D_{x-}^{\alpha} u(x, y, t) &= \displaystyle \frac{1}{\varGamma(m-\alpha)}\left(-\frac{\partial}{\partial x}\right)^{m} \int_{x}^{b} \frac{u(\xi, y, t)}{(\xi-x)^{\alpha+1-m}} d \xi \\ D_{y+}^{\alpha} u(x, y, t) &= \displaystyle \frac{1}{\varGamma(m-\alpha)}\left(\frac{\partial}{\partial y}\right)^{m} \int_{c}^{y} \frac{u(x, \eta, t)}{(y-\eta)^{\alpha+1-m}} d \eta \\D_{y-}^{\alpha} u(x, y, t) &= \displaystyle \frac{1}{\varGamma(m-\alpha)}\left(-\frac{\partial}{\partial y}\right)^{m} \int_{y}^{d} \frac{u(x, \eta, t)}{(\eta-y)^{\alpha+1-m}} d \eta \end{aligned} 其中 m-1<\alpha\leqslant mm 为正整数,c\left(u,|\cdot|\right) 为四个方向的扩散系数.
\quad\quad通过 Grünwald-Letnikov 导数对 Riemann-Liouville 导数进行离散化. 一般地,\begin{array}{l}D_{x_{i}+}^{\alpha} u\left(x_{i}, y_{j}, t_{k}\right) \approx \displaystyle \frac{1}{h^{\alpha}}\sum_{m=0}^{i+1} g_{m}^{\alpha} u\left(x_{i+1-m},y_{j},t_{k}\right),\\[1ex]D_{x_{i}-}^{\alpha}u\left(x_{i}, y_{j}, t_{k}\right) \approx \displaystyle \frac{1}{h^{\alpha}}\sum_{m=0}^{M+1-i}g_{m}^{\alpha}u\left(x_{i+m},y_{j},t_{k}\right),\\[1ex]D_{y_{i}+}^{\alpha}u\left(x_{i},y_{j},t_{k}\right) \approx \displaystyle \frac{1}{h^{\alpha}} \sum_{n=0}^{j+1} g_{n}^{\alpha}u\left(x_{i}, y_{j+1-n},t_{k}\right),\\[1ex]D_{y_{i}-}^{\alpha} u\left(x_{i}, y_{j}, t_{k}\right) \approx \displaystyle \frac{1}{h^{\alpha}}\sum_{n=0}^{N+1-j} g_{n}^{\alpha} u\left(x_{i}, y_{j+n}, t_{k}\right),\end{array} 其中 x_i=a+ih_x, (i=0,1,2,\dots,M+1)y_i=c+jh_y, (j=0,1,2,\dots,N+1),那么,空间步长为 h_x=(b-a)/(M+1)h_y=(b-a)/(N+1),时间步长为 \tau = T/Kk=1,2,\dots,K,这里 g_{m}^{\alpha}=(-1)^{m}\left(\begin{matrix}\alpha\\m\end{matrix}\right)=(-1)^{m}\frac{\alpha(\alpha-1) \cdots(\alpha-m+1)}{m !}. \quad\quad至此可以给出方程的显式有限差分格式:\begin{aligned}\frac{u_{i, j}^{k+1}-u_{i, j}^{k}}{\tau}=& \frac{\left.q_{x}\right|_{x_{i}, y_{j}, t_{k}}-\left.q_{x}\right|_{x_{i-1}, y_{j}, t_{k}}}{h}-\frac{\left.q_{y}\right|_{x_{i}, y_{j}, t_{k}}-\left.q_{y}\right|_{x_{i}, y_{j-1}, t_{k}}}{h} \\=&\frac{1}{h^{\alpha+1}}\left(c_{i, j}^{x,+} \sum_{m=1}^{i+1} g_{i+1-m}^{\alpha} u_{m, j}^{k}-c_{i-1, j}^{x,+} \sum_{m=1}^{i} g_{i-m}^{\alpha} u_{m, j}^{k}\right) \\&+\frac{1}{h^{\alpha+1}}\left(c_{i-1, j}^{x,-} \sum_{m=i-1}^{M} g_{m+1-i}^{\alpha} u_{m, j}^{k}-c_{i, j}^{x,-}\sum_{m=i}^{M} g_{m-i}^{\alpha} u_{m, j}^{k}\right) \\&+\frac{1}{h^{\alpha+1}}\left(c_{i, j}^{y,+}\sum_{n=1}^{j+1} g_{j+1-n}^{\alpha} u_{i, n}^{k}-c_{i, j-1}^{y,+} \sum_{n=1}^{j} g_{j-n}^{\alpha} u_{i,n}^{k}\right) \\&+\frac{1}{h^{\alpha+1}}\left(c_{i, j-1}^{y,-} \sum_{n=j-1}^{N} g_{n+1-j}^{\alpha} u_{i,n}^{k}-c_{i, j}^{y,-} \sum_{n=i}^{N} g_{n-j}^{\alpha} u_{i, n}^{k}\right).\end{aligned} \quad\quad接下来就尝试给出这个格式的矩阵形式,首先必须注意到,在零边值条件下,所有 u_{0,j}^{k},u_{M+1,j}^{k},u_{i,0}^{k},u_{i,N+1}^{k} 的值均为 0,因此仅需考虑 1 \leqslant i \leqslant M, 1 \leqslant j \leqslant N 时的 u_{i,j}^{k} 值. 接下来定义两个方向的离散矩阵:\mathbf{A}_x=\left(\begin{matrix} g_{1}^{\alpha} & g_{0}^{\alpha} & 0 & \cdots & \cdots & 0 \\g_{2}^{\alpha}&g_{1}^{\alpha} & g_{0}^{\alpha} & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots &\vdots \\ \vdots &\ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & g_{1}^{\alpha} &g_{0}^{\alpha} \\g_{M}^{\alpha} & \cdots & \cdots & \cdots & g_{2}^{\alpha} & g_{1}^{\alpha} \end{matrix}\right),\mathbf{B}_x=\left(\begin{matrix}g_{0}^{\alpha} & 0 & 0 & \cdots & \cdots & 0 \\ g_{1}^{\alpha} & g_{0}^{\alpha}& 0 & \ddots &\ddots& 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots &\ddots & \ddots & 0 \\ \vdots &\ddots & \ddots & \ddots&g_{0}^{\alpha} & 0 \\ g_{M-1}^{\alpha} & \cdots & \cdots & \cdots& g_{1}^{\alpha} &g_{0}^{\alpha} \end{matrix}\right), \mathbf{A}_y=\left(\begin{matrix} g_{1}^{\alpha} & g_{0}^{\alpha} & 0 & \cdots & \cdots & 0 \\g_{2}^{\alpha}&g_{1}^{\alpha} & g_{0}^{\alpha} & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots &\vdots \\ \vdots &\ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & g_{1}^{\alpha} &g_{0}^{\alpha} \\g_{N}^{\alpha} & \cdots & \cdots & \cdots & g_{2}^{\alpha} & g_{1}^{\alpha} \end{matrix}\right),\mathbf{B}_y=\left(\begin{matrix}g_{0}^{\alpha} & 0 & 0 & \cdots & \cdots & 0 \\ g_{1}^{\alpha} & g_{0}^{\alpha}& 0 & \ddots &\ddots& 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots &\ddots & \ddots & 0 \\ \vdots &\ddots & \ddots & \ddots&g_{0}^{\alpha} & 0 \\ g_{N-1}^{\alpha} & \cdots & \cdots & \cdots& g_{1}^{\alpha} &g_{0}^{\alpha} \end{matrix}\right). \quad\quad现在关注差分格式右侧的第一项,有 c_{i, j}^{x,+} \sum_{m=1}^{i+1} g_{i+1-m}^{\alpha} u_{m, j}^{k}=\left(\text{diag} \left(c_{1, j}^{x,+}, c_{2, j}^{x,+}, \cdots, c_{M, j}^{x,+}\right)\mathbf{A}_x\left(u_{1, j}^{k}, u_{2, j}^{k}, \cdots, u_{M, j}^{k}\right)^{\mathrm{T}}\right)_i. \quad\quad出于这种考虑,定义两个分块矩阵 :\mathbf{A}=\left(\mathbf{A}_1,\mathbf{A}_2,\cdots,\mathbf{A}_N \right)\mathbf{B}=\left(\begin{matrix}\mathbf{B}_1\\\mathbf{B}_2\\\cdots\\\mathbf{B}_M\end{matrix} \right),其中 \begin{aligned}\mathbf{A}_j = \ & \text{diag} \left(c_{1, j}^{x,+}, c_{2, j}^{x,+}, \cdots, c_{M, j}^{x,+}\right)\mathbf{A}_x-\text{diag} \left(c_{0, j}^{x,+}, c_{1, j}^{x,+}, \cdots, c_{M-1, j}^{x,+}\right)\mathbf{B}_x\\&+\text{diag} \left(c_{0, j}^{x,-}, c_{2, j}^{x,-}, \cdots, c_{M-1, j}^{x,-}\right)\mathbf{A}_x^{\mathrm{T}}-\text{diag} \left(c_{1, j}^{x,-}, c_{1, j}^{x,-}, \cdots, c_{M, j}^{x,-}\right)\mathbf{B}_x^{\mathrm{T}}\end{aligned} \begin{aligned}\mathbf{B}_i = \ & \text{diag} \left(c_{i,1}^{y,+}, c_{i,2}^{y,+}, \cdots, c_{i,N}^{y,+}\right)\mathbf{A}_y-\text{diag} \left(c_{i,0}^{y,+}, c_{i,1}^{y,+}, \cdots, c_{i, N-1}^{y,+}\right)\mathbf{B}_y\\&+\text{diag} \left(c_{i, 1}^{y,-}, c_{i, 2}^{y,-}, \cdots, c_{i, N-1}^{y,-}\right)\mathbf{A}_y^{\mathrm{T}}-\text{diag} \left(c_{i, 1}^{y,-}, c_{i, 2}^{y,-}, \cdots, c_{i, N}^{y,-}\right)\mathbf{B}_y^{\mathrm{T}},\end{aligned} 再对矩阵\left(u_{i,j}^k\right) 进行两种分块:\left(u_{i,j}^k\right)=\left(\mathbf{u}_{:,1}^k, \mathbf{u}_{:,2}^k, \cdots, \mathbf{u}_{:,N}^k\right)=\left(\begin{matrix}\mathbf{u}_{1,:}^k\\\mathbf{u}_{2,:}^k\\\vdots\\\mathbf{u}_{M,:}^k\end{matrix}\right),\mathbf{P}=\left(\mathbf{p}_1,\mathbf{p}_2,\cdots,\mathbf{p}_N\right),\mathbf{S}=\left(\begin{matrix}\mathbf{s}_{1}\\\mathbf{s}_{2}\\\vdots\\\mathbf{s}_{M}\end{matrix}\right),最终就有 \mathbf{p}_j=\mathbf{A}_j\mathbf{u}_{:,j}^k,\mathbf{s}_i=\mathbf{u}_{i,:}^k\mathbf{B}_i,那么 \left(u_{i,j}^{k+1}\right)=\left(u_{i,j}^{k}\right)+\tau(\mathbf{P}+\mathbf{S}). \quad\quad以上的矩阵形式已经能够满足对数值求解的要求. 但是为了使用快速显式扩散算法,需要将差分格式整理成 \mathbf{U}_{k+1}=\left(\mathbf{I}+\tau\mathbf{Q}\right)\mathbf{U}_{k} 的形式,并给出矩阵 \mathbf{Q} 谱半径的上界. 这里 \mathbf{U}_{k}\left(u_{i,j}^k\right) 的列展开,\mathbf{U}_{0}=\mathbf{F} 由初值条件给出. 现在将差分格式整理成 \begin{aligned}\frac{u_{i, j}^{k+1}-u_{i, j}^{k}}{\tau}h^{\alpha+1} =& \sum_{m = 1}^{i-2}\left(c_{i, j}^{x,+} g_{i+1-m}^{\alpha}-c_{i-1, j}^{x,+} g_{i-m}^{\alpha}\right)u_{m, j}^{k}\\&+\left(c_{i, j}^{x,+} g_{2}^{\alpha}-c_{i-1, j}^{x,+} g_{1}^{\alpha}+c_{i-1, j}^{x,-} g_{0}^{\alpha}\right)u_{i-1, j}^{k}\\&+\left(c_{i, j}^{x,+} g_{1}^{\alpha}-c_{i-1, j}^{x,+} g_{0}^{\alpha}+c_{i-1, j}^{x,-} g_{1}^{\alpha}-c_{i, j}^{x,-} g_{0}^{\alpha}\right)u_{i, j}^{k}\\&+\left(c_{i-1, j}^{x,-} g_{2}^{\alpha}-c_{i, j}^{x,-} g_{1}^{\alpha}+c_{i, j}^{x,+} g_{0}^{\alpha}\right)u_{i+1, j}^{k}\\&+\sum_{m = i+2}^{M}\left(c_{i-1, j}^{x,-} g_{m+1-i}^{\alpha}-c_{i, j}^{x,-} g_{m-i}^{\alpha}\right)u_{m, j}^{k}\\&+\sum_{n = 1}^{j-2}\left(c_{i, j}^{y,+} g_{j+1-n}^{\alpha}-c_{i-1, j}^{y,+} g_{j-n}^{\alpha}\right)u_{i, n}^{k}\\&+\left(c_{i, j}^{y,+} g_{2}^{\alpha}-c_{i, j-1}^{y,+} g_{1}^{\alpha}+c_{i, j-1}^{y,-} g_{0}^{\alpha}\right)u_{i, j-1}^{k}\\&+\left(c_{i, j}^{y,+} g_{1}^{\alpha}-c_{i, j-1}^{y,+} g_{0}^{\alpha}+c_{i, j-1}^{y,-} g_{1}^{\alpha}-c_{i, j}^{y,-} g_{0}^{\alpha}\right)u_{i, j}^{k}\\&+\left(c_{i, j-1}^{y,-} g_{2}^{\alpha}-c_{i, j}^{y,-} g_{1}^{\alpha}+c_{i, j}^{y,+} g_{0}^{\alpha}\right)u_{i, j+1}^{k}\\&+\sum_{n = j+2}^{N}\left(c_{i, j-1}^{y,-} g_{n+1-j}^{\alpha}-c_{i, j}^{y,-} g_{n-j}^{\alpha}\right)u_{i, n}^{k},\end{aligned} 注意上式右侧的后两项,考虑两个分块 N \times N 矩阵 \mathbf{Q}_{x}\mathbf{Q}_{y} ,它们每一块形状均为 M \times M ,其中 \mathbf{Q}_{x} 是一个分块对角矩阵,其中每一块 \left(\mathbf{Q}_{x}\right)_j 均有如下形式:

\frac{1}{h^{\alpha+1 }} \left(\begin{matrix}c_{1, j}^{x,+} g_{1}^{\alpha}-c_{0, j}^{x,+} g_{0}^{\alpha}+c_{0, j}^{x,-} g_{1}^{\alpha}-c_{1, j}^{x,-} g_{0}^{\alpha}&c_{0, j}^{x,-} g_{2}^{\alpha}-c_{1, j}^{x,-}g_{1}^{\alpha} +c_{1, j}^{x,+} g_{0}^{\alpha}&c_{0, j}^{x,-} g_{3}^{\alpha}-c_{1, j}^{x,-}g_{2}^{\alpha}& \cdots & c_{0, j}^{x,-} g_{M}^{\alpha}-c_{1, j}^{x,-}g_{M-1}^{\alpha}\\c_{2, j}^{x,+} g_{2}^{\alpha}-c_{1, j}^{x,+}g_{1}^{\alpha} +c_{1, j}^{x,-} g_{0}^{\alpha}&c_{2, j}^{x,+} g_{1}^{\alpha}-c_{1, j}^{x,+} g_{0}^{\alpha}+c_{2, j}^{x,-} g_{1}^{\alpha}-c_{1, j}^{x,-} g_{0}^{\alpha}&c_{1, j}^{x,-} g_{2}^{\alpha}-c_{2, j}^{x,-}g_{1}^{\alpha} +c_{2, j}^{x,+} g_{0}^{\alpha}& \cdots & c_{1, j}^{x,-} g_{M-1}^{\alpha}-c_{2, j}^{x,-}g_{M-2}^{\alpha}\\c_{3, j}^{x,+} g_{3}^{\alpha}-c_{2, j}^{x,+}g_{2}^{\alpha}&c_{3, j}^{x,+} g_{2}^{\alpha}-c_{2, j}^{x,+}g_{1}^{\alpha} +c_{2, j}^{x,-} g_{0}^{\alpha}&c_{3, j}^{x,+} g_{1}^{\alpha}-c_{2, j}^{x,+} g_{0}^{\alpha}+c_{3, j}^{x,-} g_{1}^{\alpha}-c_{2, j}^{x,-} g_{0}^{\alpha}& \cdots & c_{2, j}^{x,-} g_{M-2}^{\alpha}-c_{3, j}^{x,-}g_{M-3}^{\alpha}\\\vdots & \ddots & \ddots & \ddots &\vdots\\c_{M, j}^{x,+} g_{M}^{\alpha}-c_{M-1, j}^{x,+}g_{M-1}^{\alpha}& \cdots & \cdots & \cdots &c_{M, j}^{x,+} g_{1}^{\alpha}-c_{M-1, j}^{x,+} g_{0}^{\alpha}+c_{M-1, j}^{x,-} g_{1}^{\alpha}-c_{M, j}^{x,-} g_{0}^{\alpha}\end{matrix}\right),

\mathbf{Q}_{y} 的每一块都是一个对角矩阵,其表达较上式更加繁琐,在这里不再给出其完整的矩阵表示,而是直接给出 \mathbf{Q}=\mathbf{Q}_{x}+\mathbf{Q}_{y} 的定义: \begin{aligned}\left(Q_{j, j}\right)_{i, m}=& \frac{1}{h^{\alpha+1}}\left(c_{i, j}^{x,+} g_{i+1-m}^{\alpha}-c_{i-1, j}^{x,+} g_{i-m}^{\alpha}\right), \quad m \leqslant i-2 \\\left(Q_{j, j}\right)_{i, i-1}=& \frac{1}{h^{\alpha+1}}\left(c_{i, j}^{x,+} g_{2}^{\alpha}-c_{i-1, j}^{x,+} g_{1}^{\alpha}+c_{i-1, j}^{x,-} g_{0}^{\alpha}\right) \\\left(Q_{j, j}\right)_{i, i}=& \frac{1}{h^{\alpha+1}}\left(c_{i, j}^{x,+} g_{1}^{\alpha}-c_{i-1, j}^{x,+} g_{0}^{\alpha}+c_{i-1, j}^{x,-} g_{1}^{\alpha}-c_{i, j}^{x,-} g_{0}^{\alpha}\right.\\&\left.+c_{i, j}^{y,+} g_{1}^{\alpha}-c_{i, j-1}^{y,+} g_{0}^{\alpha}+c_{i, j-1}^{y,-} g_{1}^{\alpha}-c_{i, j}^{y,-} g_{0}^{\alpha}\right) \\\left(Q_{j, j}\right)_{i, i+1}=& \frac{1}{h^{\alpha+1}}\left(c_{i-1, j}^{x,-} g_{2}^{\alpha}-c_{i, j}^{x,-} g_{1}^{\alpha} + c_{i, j}^{x,+} g_{0}^{\alpha}\right) \\\left(Q_{j, j}\right)_{i, m}=& \frac{1}{h^{\alpha+1}}\left(c_{i-1, j}^{x,-} g_{m+1-i}^{\alpha}-c_{i, j}^{x,-} g_{m-i}^{\alpha}\right), \quad m \geqslant i+2 \\\left(Q_{j, n}\right)_{i, m}=& \frac{1}{h^{\alpha+1}} \delta_{i, m}\left(c_{i, j}^{y,+} g_{j+1-n}^{\alpha}-c_{i, j-1}^{y,+} g_{j-n}^{\alpha}\right), \quad n \leqslant j-2 \\\left(Q_{j, j-1}\right)_{i, m}=& \frac{1}{h^{\alpha+1}} \delta_{i, m}\left(c_{i, j}^{y,+} g_{2}^{\alpha}-c_{i, j-1}^{y,+} g_{1}^{\alpha}+c_{i, j-1}^{y,-} g_{0}^{\alpha}\right) \\\left(Q_{j, j+1}\right)_{i, m}=& \frac{1}{h^{\alpha+1}} \delta_{i, m}\left(c_{i, j-1}^{y-} g_{2}^{\alpha}-c_{i, j}^{y,-} g_{1}^{\alpha}+c_{i, j}^{y,+} g_{0}^{\alpha}\right) \\\left(Q_{j, n}\right)_{i, m}=& \frac{1}{h^{\alpha+1}} \delta_{i, m}\left(c_{i, j-1}^{y,-} g_{n+1-j}^{\alpha}-c_{i, j}^{y,-} g_{n-j}^{\alpha}\right), \quad n \geqslant j+2.\end{aligned} \quad\quad其中 \delta_{i, m} 为 Kronecker 符号. 至此已经完成了方程离散化的全部工作.

发表评论

您的电子邮箱地址不会被公开。