求解分数阶扩散方程的循环预处理的极小化残量法

求解分数阶扩散方程的循环预处理的极小化残量法
屈威
【摘要】将循环预处理的极小化残量法应用到分数阶扩散方程的求解中,利用Crank-Nicolson方法给出了扩散方程的隐差分格式,以及循环预处理矩阵的形式,
并通过数值实验说明.
【期刊名称】《韶关学院学报》
【年(卷),期】2016(037)004
【总页数】5页(P1-5)
【关键词】分数阶扩散方程;Toeplitz矩阵;循环预处理算子;极小残量法
【作者】屈威
【作者单位】韶关学院数学与统计学院,广东韶关512005
【正文语种】中文
【中图分类】O241.82
分数阶微分方程是传统整数阶微分方程的推广，其产生于一些反常扩散模型.分数
阶微分方程主要包括空间分数阶微分方程、时间分数阶微分方程、空间-时间分数
阶微分方程3大类.由于分数阶微分算子与整数阶微分算子相比，其具有非局部的
性质，非常适合描述现实世界中具有记忆以及遗传性质的材料，现已成为描述各类力学与物理行为的重要工具，因此被广泛应用于流体力学、模拟地下水的传送、反常扩散、信号处理与系统识别等其他科学研究领域.然而，分数阶微分方程的解析
解很难显示表出，因而分数阶微分方程的数值解与数值算法备受研究专家者的关注，
现已成为当前研究领域的热门问题.
M.M.Meerschaert等首先提出了利用有限差分的方法求解分数阶微分方程［1-2］.近年来，有关分数阶微分方程的数值解的研究工作也取得了许多显著的成果，主要包括：（1）有限差分法；（2）有限元法；（3）预处理法；（4）多重网格法；（5）反积分法等［2-9］.本文主要研究空间分数阶扩散微分方程的数值方法，运用预处理方法及极小残量法对该问题进行求解.
本文考虑如下带有Dirichlet齐次边界条件的分数阶微分扩散方程：
其中d+，d-为非负常数，且d++d-≠0，f（x，t）为定义在［a，b］×［0，T］上的已知函数，微分算子）和（x，t）分别为函数u（x，t）的α（1＜α＜2）阶
的左Riemann-Liouville和右Riemann-Liouville型分数阶导数，在文献［9］中，给出了左、右Riemann-Liouville型分数阶导数的一般定义.
设｛u（xi，tm）|0≤i≤N+1，0≤m≤M｝为微分方程（1）的解析解，
｛uim|0≤i≤N+1，0≤m≤M｝为分数阶微分方程（1）的数值解.
令uim≈u（xi，tm），fim+1/2=f（xi，tm+1/2），，利用Crank-Nicolson隐差分格式和文献［9］给出2阶左、右Riemann-Liouville的分数阶导数的微分算子，得到如下差分格式：
该差分格式是无条件稳定的，以及在空间和时间上的误差都满足2阶精度要求［10］.令:
则差分格式（2）的矩阵形式表示为：
其中，I为单位矩阵，Qa和QTa为Toeplitz
矩阵，
定义1一个N×N阶矩阵T满足:
其中tij=ti-j，这样的矩阵称为Toeplitz矩阵，显然Toeplitz矩阵的主对角线上各元素彼此相等，且平行于主对角线的各对角线上的元素也彼此相等［11］.
序列｝的表达式满足如下关系和性质：
由于矩阵M（m+1）是非对称的Toeplitz矩阵，而极小残量法（MINRES）只适用于求解对称的方程，因此，选择利用置换矩阵将方程（3）转变为对称的方程，进而求解方程（3）的等价方程：
由此可以利用MINRES方法求解方程（4）.然而，在方程（4）中，系数矩阵YM （m+1）条件数可能比较大，这将会造成使用MINRES方法时迭代收敛过慢以及存储计算量大等缺点.在这种情况下，可利用预处理方法克服以上缺点.在文献［5］中，已经详细地研究了如何利用如下Strang型循环预处理方法对方程（3）进行
数值求解，预处理矩阵为，其中循环矩阵S（Qa）、S（QTa）的第一列分别为：通过循环预处理的方程通常可以借助快速Fourier变换进行求解，这将大大提高方程求解的效率，由文献［11］可以知道，循环矩阵C可被Fourier矩阵对角化，
即C=F*ΛCF，i为虚数单位，矩阵ΛC为对角矩阵，其对角线元素分别为循环矩阵C特征值.由于S（m+1）为正规矩阵，利用文献［1］中的思想，可得如下循环预处理矩阵：
其中为对角矩阵，对角元素为S（m+1）特征值的绝对值，即，并且和分别为由
循环矩阵S（ Q a"和S（ QT a"特征值的模组成的对角矩阵.因此，预处理条件下的同解方程为：
本文主要利用Strang型循环预处理的极小残量法（P1MINRES）、Tony Chan
型循环预处理的极小残量法（P2MINRES）求解方程（5）以及极小残量法（MINRES）求解方程（4），并对以上3种方法进行比较，见表1，表2.迭代停
止条件为其中，rk表示迭代第k次的残差，限制最大的迭代次数MAXIT=2 000，在表1和表2中用“N”表示空间节点的个数；“M”代表时间的步长数；“Error”表示方程（1）解析解与数值解的误差；“Rate”表示误差精度；“CPU（s）”表示迭代收敛所需要的时间；“Iter”代表求解方程（3）达到迭代
停止条件时所需的平均迭代的次数.
考虑分数阶扩散方程（1）：取，源项：
其边界条件为和初始条件方程（1）的精确解为:
从表1可以看出，利用循环预处理的极小残量法求解例题得到的数值解满足理论的误差与2阶精度要求，这与理论结果相符，说明的方法是正确的.从表2可以看出，求解方程（5）时，经过预处理的2种方法比没有预处理的方法在计算时间上更快，所需要的迭代次数也更少，这也更能说明该方法的有效性与实用性.除此之外，P1MINRES和P2MINRES方法能在理论上得到矩阵|S（m+1）|-1YM
（m+1）的特征值分布区间，进一步可得循环预处理的极小残量法的收敛性［1，5］.
【相关文献】
［1］Pestana J，Wathen A J.A preconditioned MINRES method for nonsymmetric Toeplitz matrices，SIAM［J］.J Matrix Anal&Appl，2015，36 （1）:273-288.
［2］Meerschaert M M，Scheffler H P，Tadjeran C.Finite difference methods for two-dimensional fractional dispersion equation［J］.J Comput Phys，2006（21）:249-261. ［3］Meerschaert M M，Tadjeran C.Finite difference approximations for fractional advection–dispersion flow equations［J］.J Comput Appl Math，2004（172）:65-77. ［4］Meerschaert M M，Tadjeran C.Finite difference approximations for two-sided space-actional partial differential equations［J］.Appl Numer Math，2006（56）:80-90.
［5］Lei S L，Sun H W.A circulant preconditioner for fractional diffusion equations［J］.J Comput Phys，2013（242）:715-725.
［6］Pang H，Sun H.Multigrid method for fractional diffusion equations［J］.J Comput Phys，2012，231（2）:693-703.
［7］Pang H，Sun H.Fast numerical contour integral method for fractional diffusion equations［J］.J of Sci Comput，2015（19）:1-26.
［8］Qu W，Lei S L，Vong S W.Circulant and skew-circulant splitting iteration for fractional advection-diffusion equations［J］.Inter J of Comput Math，2014（91）:2232-2242.
［9］Sousa E，Li C.A weighted finite difference method for the fractional diffusion
equation based on the Riemann–Liouville derivative ［J］.Appl Numer Math，2015（90）:22-37.
［10］Deng W，Chen M.Efficient numerical algorithms for three-dimensional fractional partial differential equations［J］.Journal of Computational Mathematics，2014，32（4）:371-391.
［11］Chan R H，Jin X Q.An Introduction to Iterative Toeplitz Solvers
［M］.Philadelphia:Society for Industrial and Applied Mathematics，2007.。