2014年5月 May，2014 计算数学 MATHEMATICA NUMERICA SINICA 第36卷第2期 V_0136No2 空间一时间分数阶变系数对流扩散方程 微分阶数的数值反演 ) 贾现正 张大利 李功胜 (山东理工大学理学院应用数学研究所，山东淄博255049) 池光胜 (山东凯文科技职业学院本科教育学院，济南250200) 李慧玲 (山东理工大学理学院应用数学研究所，山东淄博255049) 摘 要 考虑终值数据条件下一维空间一时间分数阶变系数对流扩散方程中同时确定空间微分阶数 与时间微分阶数的反问题基于对空间一时间分数阶导数的离散，给出求解正问题的一个隐式差 分格式，通过对系数矩阵谱半径的估计，证明差分格式的无条件稳定性和收敛性联合最佳摄动量 算法和同伦方法引入同伦正则化算法，应用一种单调下降的Sigmoid型传输函数作为同伦参数， 对所提微分阶数反问题进行精确数据与扰动数据情形下的数值反演结果表明同伦正则化算法对 于空间一时间分数阶反常扩散的参数反演问题是有效的 关键词：空间一时间分数阶扩散；差分格式；稳定性与收敛性；微分阶数反问题；同伦正则化 算法；数值反演 MR(2000)主题分类：35R30，35Rll，65M06 1引 言 近二十年来，反常扩散行为引起人们越来越多的关注，分数阶扩散方程在物理学、力学、 环境科学、水文地质学及金融学等领域得到了广泛的应用1 vJ譬如，污染物在多孔介质中的 传输往往存在着非对称的拖尾(heavy tail)或提前穿透(early breakthrough)现象，这就是溶 质运移的非费克(NonFickian)反常扩散对于这类反常扩散行为，已有研究表明，分数阶扩 散模型相比经典的整数阶高斯(Gauss)模型能够更好地模拟扩散过程和重建试验数据分数 阶导数具有记忆性、遗传性和整体性，特别对于污染物长距离传输时表现出的非对称、非线 性的整体性行为模式的研究，分数阶扩散模型可能是一个很有效的研究工具通常，若考虑时 间相关性或扩散的记忆性，就得到时间分数阶扩散方程，其中时间微分阶数介于0与1之间 的反常扩散要比正常扩散慢；若考虑空间相关性或非局域性，则得到空间分数阶扩散方程，其 中空间微分阶数介于1与2之间的反常扩散相比正常扩散要快一些如果既考虑时间相关性， 又考虑空间整体性，则得到所谓的空间一时间分数阶扩散方程 2013年5月21日收到 )基金项目：国家自然科学基金资助项目(11071148)和山东省自然科学基金资助项目(ZR2011AQ014) 。)通讯作者：李功胜，Email：ligssduteducn 114 计算数学 2014钲 上世纪90年代以来，相继出版的专著8-1oj对于分数阶导数的概念，若干分数阶微分方 程解的适定性理论分析、解析解与数值解法，以及分数阶微分方程在科学与工程中的应用研 究等给予了综述总结，这其中关于时间分数阶扩散的研究居多在时间分数阶扩散方程的应 用研究方面，Adams与Gelhar【 以及Hatano等L2_较早发表了有影响的研究成果，揭示了实 际存在的一些反常扩散行为在时间分数阶微分方程初边值问题理论分析研究方面，Gorenflo 等11_ ，Hanygad13，Eidelman与Kochubei14，以及Luchko15 17】等作出了开创性的研 究，证明了解的存在唯一性，以及关于时间分数阶扩散的极值原理等结论对于时间分数阶 扩散方程的数值方法，刘发旺(FLiu)及其研究团队18-2oJ，Murio ，国内许传炬团队_22_， 以及孙志忠团队l 2 3_等主要在差分法研究方面做出了贡献对于空间分数阶扩散方程的研究， Meerschaert和Benson研究团队的工作24-28 J最具影响，刘发旺及其研究团队的工作29 32 J 也广受关注，国内陈文等 3 3_、邓伟华等_34_35在数值解法及应用方面也有较好的研究 另一方面，一旦建立了一个分数阶扩散模型，求解和分析模型中难以直接测量的关键参 数就成为一类很重要的科学问题比如，反映空间或时间相关性的微分阶数的确定问题，反映 空间非均质性的扩散系数的确定问题，反映外界力作用属性的源汇项的确定问题，等等这 就导致了对于分数阶反常扩散中的反问题研究虽然对于分数阶扩散中的反问题研究还不多 见，但近年来一些工作已引起了大家的关注这其中对于时间分数阶扩散方程反问题的研究 较多，见Murio36J，程晋等 ，Bondarenko与Ivaschenkoa8 J，刘继军和Yamamoto39J，郑 光辉与魏婷_4o_，Sakamoto与Yamamoto41，徐翔等 ，Tuan4a】，Yamamoto等 ，Jin与 Rundell45】，李功胜等【46_ ，熊向团等【48，魏婷等49此外，对于空间分数阶扩散模型，由于 正问题研究的困难性，相应反问题的研究相对较少，见韦慧等【50，郑光辉与魏婷 ，池光胜 等_I52，邓伟华等 随着对于分数阶扩散正演反演问题的研究进展，有两个方面值得关注：一是分数阶扩散 在空间多维数情形下的相关问题研究目前对于分数阶高维扩散，反问题研究很少；而鉴于空 间分数阶导数的复杂性和计算复杂度，高维正问题研究也不是很多第二个方面是对于既含 时间分数阶导数、又含空间分数阶导数的空间一时间分数阶扩散相关问题的研究，这类反常 扩散模型的研究要比单纯的时间分数阶或空间分数阶扩散更难对于(左侧)空间一时间分数 阶扩散正问题的求解，文19，54，55等曾给出了差分求解格式，并证明了格式的稳定性和收 敛性 本文继续考虑(左侧)空间一时问分数阶变系数对流扩散方程正问题的差分解法，同时将 探讨根据终值数据确定两个微分阶数的反演问题对于正问题的差分格式，文中将通过对系 数矩阵谱半径的估计，证明格式的无条件稳定性和收敛性这种证明方法相比已有方法显得更 加简洁清晰，且可应用于一般时间空间分数阶扩散方程差分格式的理论分析注意到微分阶 数是刻划分数阶扩散的首要指标，也是分数阶扩散区别于整数阶Gauss扩散的重要参数由 于实际问题中微分阶数是未知的，更是难以通过实验手段直接测量获得，因而对于分数阶扩散 中微分阶数的反问题研究具有重要意义目前，除了程晋等_37J，Bondarenko与Ivaschenko38J 对于一维时间分数阶扩散中微分阶数与扩散系数联合反演问题的研究之外，尚未见有其他文 献报道关于微分阶数反问题的研究 本文将从数值方法角度，探讨一维(左侧)空间一时间反常扩散模型中两个微分阶数的数 值反演问题文中将采用一种结合同伦思想与最佳摄动量算法的同伦正则化算法从已有文 献看，最佳摄动量算法无论对于整数阶扩散模型的参数反演，还是对于分数阶反常扩散模型的 参数反演问题求解都取得了很好的数值结果，见文【46，47，50，52，5659等然而，对于同时确 2期 贾现正等：空间一时间分数阶变系数对流扩散方程微分阶数的数值反演 115 定多个不同类型参数的病态性更强的多参数反演问题，常规的最佳摄动量算法可能会失去效 用，我们需要更精细地选取正则参数以克服反演的严重病态性另一方面，同伦方法是研究非 线性代数方程求解的有力工具，而联合同伦思想与摄动方法的所谓同伦摄动方法(homotopy perturbation method)在非线性微分方程求解研究方面发挥了重要作用，见文6062】等对于 分数阶反常扩散方程的多参数反演，一种可行的途径是联合同伦方法与最佳摄动量算法，实施 同伦正则化算法这种迭代算法可以扩大收敛区域，且已被应用于整数阶扩散相关反问题的 数值求解，见【63，64等不过，本文仍然应用这种同伦正则化算法，主要基于以下两个原因： (i)正问题更为困难文中要求解空间时间分数阶反常扩散方程，计算复杂度较高； (ii)反演问题不同文中要考虑两个微分阶数的同时反演问题，这要比先前的工作病态性 更高 本文余下的主要内容安排如下 第2节给出求解一维空间时间分数阶变系数对流扩散方程正问题的隐式差分格式，并 证明格式的稳定性与收敛性；第3节对于由终值数据确定微分阶数的反问题，讨论一种以Sig moid型传输函数为同伦参数的同伦正则化反演算法；第4节分别在精确数据与扰动数据情形 下，应用同伦正则化算法对微分阶数进行数值反演；最后在第5节进行总结，给出本文的主要 结论 2正问题及其数值求解 对于f0，T0，考虑一维空间一时间分数阶变系数对流扩散方程 一 a ( ) +s( ，)，。0为 扩散系数，v(x)0为对流系数，s( ，t)为源汇项显见，若OL=1， =2，则方程(1)为通常 整数阶的经典对流扩散方程方程(1)中的时间分数阶导数用Caputo的定义： ： 1 一 at r( 一 ) 0 。 a丁 、 而空间分数阶导数用RiemannLiouville的定义： = r(2 Ox2 武 az1 一吖) 。n (z一)11一 。 、 对于方程(1)，给定初始条件为 u(x，0)=u0( )，0 z， 边值条件为 u(o，t)=u(1，t)=0，0 0， ：0，1，M；t ： 丁，丁0， =0，1，其中h和7-分别是空间和时间步长记 札 =u( t，t )，s =s( i，tk)，在(z ，tk+1)处，Caputo分数阶导数离散为 0 u( ，tk+1) Ota (u “一l u 一。)【 +1) 一。一j1-+0(下) =F (，“ + 一utk一乱 【( +1) 一 一kl- + f“ kj(J+2) 一。+J 一a一2(j+1) 一a)+O(T) +f“t一 +2) 一。+ 一。一 +1) 一。)+ ) J=1 而Riemann-Liouville分数阶导数采用移位的Grfinwald公式离散为 Ou(xi，tk+1) 一= a 方程中的一阶对流项用通常的一阶向后差分近似为 ou( ，tk+1) 札 k+ 一 uu k一+1l Ox h +O(h) 记 “ =“o(zi)，s =s(xi，tk)，以及Di=J)( i)，Vi=(z )，并记 f ：2j 一 一 +1) 一 一 1) 一 ，J=1，2，尼， l bk：( +1) 一 一 一 ， =1，2，一1 注意到边界条件 6=u =0，方程(1)离散为 再记 k+1 +1 =一 整理可得 当 一0时 Vi 一1 一bku+dju -j) J=l + 1 i+ Di 7- r(2 u 1+(1+Pirig1)U 当 0时 】 一，fpt u +1+stk+l，i=1，2，M一1 ， h )JJ r(j一 ) r(-7)r(j+1) i+1 ( +哪2) 一1一 一=u 0+7- r(2一 )s 1 j=3 k +l(1 钆 一(Pi+rig2 k+l 薹 “ ：u (22 一a)+ + 一j+b 钆 +7-nr(2一a)s +1， j=1 (8) 2期 贾现正等：空间一时间分数阶变系数对流扩散方程微分阶数的数值反演 1 17 ： 薹 (14) (15) 22差分格式的收敛性与稳定性 本小节证明上述差分格式(14)的系数矩阵 严格对角占优，且差分格式是无条件稳定 和收敛的为此，需要下述两个引理 引理1若7(1，2)，则对于(11)式中给出的 成立 gj=0； (16) 其中go=1，gl=一，y，且对于所有J 2，gj0 证明由(11)式，改写夕J为gj：(一1) ( )，对任意70，注意到下述展开式 (1 k ( J 【0 J 2 引理2对于(9)式定义的 6 ，成立 dj+bk=1， =1 2一，一1 证明根据(9)式，对于k=1，2，一1，