两项线性分数阶微分方程的奇点分离分段配置法

李玉玉，廉 欢，王同科

（天津师范大学数学科学学院，天津 300387）

分数阶微分方程是分数阶微积分的重要分支，其中，含有Caputo 分数阶导数的微分方程可以有效地描述数学、物理及工程等领域的一些实际问题.一般情况下，该类微分方程的精确解很难解析求得，故而利用数值方法求其高精度的近似解具有重要的研究价值.常用的数值方法有Laplace 变换法[1]、有限差分方法[2]、外推法[3]和Adomian 分解法[4]等.

考虑具有两项Caputo 型导数的线性分数阶微分方程

其中：0 ＜β ＜α，1 ＜α≤2；c1、c2为常数；强制项f（x）为在零点具有奇性的已知函数，且其在x=0 处存在Puiseux 级数展开式

其中[β]为β 的整数部分.文献[10]对于α/β 为有理数的情形，给出了方程（1）（仅对齐次方程f（x）=0）的解稳定的充分必要条件.文献[11]给出了阶数均小于等于2 的三项齐次方程的解稳定的充分必要条件.方程（1）可以使用Laplace 变换法求解，但求其逆变换时将遇到很大困难，文献[12]提出逆分数阶Shehu 变换方法来求其精确解，文献[13-14]设计了逆Laplace 变换法求得相应方程的渐近解.求解方程（1）的另一类思路[15]是利用其等价形式（4），积分方程（4）具有弱奇异核，其解u（x）通常在x=0 处二阶导数奇异，导致传统数值方法的整体计算精度降低，例如，直接在全区间利用分段插值函数做配置法，其收敛阶与解的奇性保持一致[16].为进一步提高计算精度，文献[16-17]采用渐变网格求解标准的Volterra 弱奇异积分方程，得到了在分段多项式空间的配置解，并给出了误差阶估计.对于渐变网格，算法的阶数与渐变指数有关，当网格在零点附近过密时，由于舍入误差的影响往往很难达到理论的收敛阶.文献[18]利用函数在奇点的Puiseux级数，通过Picard 迭代与符号运算得出第二类弱奇异Volterra 积分方程的解在零点的渐近展开式，准确地刻画了解的奇性，从而在较大区间内可以利用数值方法求解方程，有效降低了奇性对后续区间数值收敛阶的影响.

本文设计一种新的数值方法求解方程（1），对于其等价形式（4），首先使用Picard 迭代求出方程的解在零点的有限项Puiseux 级数展开式，在包含零点的小区间上利用该级数展开式代替解进行运算，将奇异方程转化为正则方程；然后在较大区间上利用Lagrange分段插值求方程的配置解，并进行收敛性分析.最后，通过数值算例验证所提方法的有效性.

1 解在零点的渐近展开式

由此可得如下定理1.

定理1 设方程（1）的强制项f（x）在x=0 处存在级数展开式（2），通过Picard 迭代求解相应的积分方程（5），则迭代收敛到解y（x）在x=0 处的Puiseux 级数展开式y∞（x），即方程（1）的解u（x）在x=0 处有Puiseux级数展开式u∞（x）=u0+y∞（x），其中：βi*为实数且满足1≤β1*＜β2*＜…→+∞，pi*、v*i，j均为非负整数.

在实际计算中，由式（9）可知，随着Picard 迭代的进行，新增项关于x 的幂次不断增大.通过控制迭代次数ρ，可保证yρ（x）中有足够多的项保持不变，进而可得到方程（5）的解y（x）在x=0 附近具有较高精度的有限项级数展开式. 基于此，针对方程（4）设计如下Puiseux 级数展开算法.

第1 步：选取适当大的正整数M，表示将要恢复出的级数展开式中x 的最高幂次.令y0（x）=g（x）满足式（8），因取有限项，要求αi+α≤M.

第2 步：由前面的讨论知，y0（x）中的xβ1，0项在以后的迭代中保持不变，对此项进行一次Picard 迭代，得

显然，yρ（x）中幂次低于β*ρ，ρ-1+α-β 的项不再改变.

第4 步：继续执行式（12）的迭代，直到β*ρ，ρ-1+αβ≥M 时停止.此时，yρ（x）中幂次低于M 的项不再改变，可得方程（5）的解y（x）在零点的渐近展开式的有限项截断为

本节通过Picard 迭代及Puiseux 级数展开式的符号运算得到了u（x）在x=0 处的有限项渐近展开式up（x）.该级数展开式具有局部逼近性质，根据定理2 可知，当x →0 时，渐近解up（x）有很高的逼近精度，但随着x 不断增大，其逼近精度将逐渐降低，因此在较大的区间上仍需要设计有效的数值方法来进一步求解方程（4）.

2 分段配置法

在方程（4）中令γ1= α - β，γ2= α，由于积分核（x-t）γσ-1（σ=1、2）中存在奇性，方程的解u（x）在x=0处二阶导数奇异.在全区间上，无论采用何种数值方法，均会因解在初始点的奇性导致全局收敛阶降低.前面导出了解u（x）在x=0 处的渐近展开式up（x），其在x=0 附近有很高的逼近精度.因此，可以分离出解的奇异部分，在剩余区间上采用分段配置法计算数值解，以提高后续区间上数值方法的计算精度.

假设方程（5）有唯一解y（x）[17，19]，下面的目标是得到区间[0，T]上的数值解，其中T 是区间上界.由于y（x）在x=0 处仅一阶可导，给定一个适当大的正数δ，将区间[0，T]划分为[0，δ]∪[δ，T]. 在[0，δ]上，使用截断Puiseux 级数展开式yp（x）来近似代替方程（5）中的y（x），以保证小区间上具有较高的计算精度. 在区间[δ，T]上采用分段配置法求解方程（5）.给定正整数N，将区间[δ，T]划分为均匀网格Ih={δ=x0＜x1＜…＜xN=T}，网格步长为h =（T - δ）/N. 在小区间ηn= [xn，xn+1]（0≤n≤N - 1）上插入m 个配置点，形成配置点集合Xh={xn，i=xn+dih：0 ＜d1＜…＜dm=1}，其中di是相应的配置参数. 对任意函数y（x）∈C[δ，T]，在Ih上构造分段m 次插值多项式πmy（x），其在ηn上的局部表达式为由于yp（x）在x=0 附近逼近y（x）的精度很高，只要选取适当的δ，可使上述逼近达到需要的精度.对于方程（17），由于y（x）在区间[δ，T]上充分光滑，任何标准数值方法均可使用，此处采用分段插值方法.

在方程（17）中使用配置解yh（x）∈Sm(0)（Ih）代替y（x），并令其在配置集合Xh上成立，则相应的配置方程为

其中yh（δ）=g～δ（δ）.

由于f（t）在t=0 处可能奇异，I0x（α；f）是一个两端点均奇异的积分，可以使用改进的复合Gauss-Legendre（GL）求积公式[20]高效计算.积分I0δ（γσ；yp）（x ＞δ）虽然仅在t=0 处有奇性，但当x 接近δ 时，被积函数在t=δ处接近奇异，也需要特殊处理，处理方法如下：当xδ≥0.2 时，yp（t）使用Padé 逼近yr（t）来近似表示，积分I0δ（γσ；yr）利用改进的GL 求积公式计算；当x-δ ＜0.2时，使用解析方法计算I0δ（γσ；yp），代入yp（x）的表达式（13），得

上式中的不完全Beta 函数的偏导数B(k，0)（δ/x，β*i+1，γσ）使用渐近展开方法来计算[21].同样，对于式（22）中的积分，当m 较小时，如m=1、2，可用解析方法直接计算，此处不再具体给出，当m 较大时，则仍需用GL 求积公式计算.

3 误差分析

当x≥δ 时，令eh=u-uh，由式（17）和式（19）可知eh满足如下误差方程：

将式（31）和式（37）代入式（30），注意到式（28），可得相应配置解的误差收敛阶，即误差估计式（29）成立.定理证毕.

由式（29）可知，在区间[δ，T]上，配置解的误差包含2 部分：在区间[0，δ]上Puiseux 级数的逼近误差和分段插值误差.显然，只要选取适当的δ，可保证配置解具有最优收敛阶.根据文献[16]提出的误差估计方法，当α-β∈（0，1）时，理论上方程（5）在全区间上的收敛阶为O（h1+α-β），在配置节点处的收敛阶可达到O（h1+2(α-β)）.本文格式的收敛阶在全区间上都可以达到O（hm+1），因此提高了计算精度.

4 数值算例

例1 考虑两项线性Caputo 分数阶微分方程初值问题

由图1 可见，渐近解up（x）、ur（x）在x = 0 处附近都有很高的逼近精度，且ur（x）在较大的区间上逼近效果更好.由于当x ≤0.6 时，up（x）仍具有O（10-14）的精度，故选取δ= 0.6，T= 20，在剩余区间[δ，T]上分别使用分段线性（m =1）和二次配置法（m =2）求所选节点的数值解un，i（i = 1，2，…，m），进一步求得相应的误差及收敛阶，分别记为E∞，δ（h）和Oδ= log2（E∞，δ（h）/E∞，δ（h/2）），并与全区间上标准的分段配置法（δ=0）的结果E∞，0（h）和O0进行比较，结果分别见表1 和表2.

图1 例1 的零点渐近解的绝对误差常用对数图像Fig.1 Absolute errors of asymptotic solutions with common logarithmic scale in Example 1

由表1 和表2 可知，渐近解up（x）和区间[δ，T]上的数值解均有较高的精度，验证了方法的有效性.由于方程的解一阶可导，[δ，T]和全区间上的线性插值格式均达到了2 阶精度；而二次插值格式在[δ，T]上达到了4 阶精度，在[0，T]上达到了2.5 阶精度，这说明分段二次格式在节点还具有超收敛性，有待进一步研究.

表1 例1 在不同步长下线性配置格式的误差及数值收敛阶Tab. 1 Errors and numerical convergence orders of linear collocationschemeunderdifferentstepsinExample1

表2 例1 在不同步长下二次配置格式的误差及数值收敛阶Tab.2 Errors and numerical convergence orders of quadratic collocation scheme under different steps in Example 1

例2 考虑两项线性Caputo 分数阶微分方程初值问题

ln x ln（1+x）/x2/3=

5.055 37×10-4x101/3-7.168 05×10-6x34+

5.109 17 ×10-6x103/3-…+6.302 95×10-6x35+

（x1/3-0.5x4/3+0.333 333x7/3-0.25x10/3+…-

1.070 88×10-4x35）ln x-ln x ln（1+x）/x2/3er（x）的具体形式不再给出，它们的绝对误差常用对数图像见图2.

图2 例2 的零点渐近解的绝对误差常用对数图像Fig.2 Absolute errors of asymptotic solutions with common logarithmic scale in Example 2

由图2 可见，渐近解up（x）和ur（x）在x=0 处附近有很高的逼近精度，但随着x 的增大，逼近精度逐渐降低.由于当x≤0.4 时，up（x）具有O（10-14）的精度，选取δ =0.4 即可保证剩余区间[δ，T]上分段配置法有较高精度的初始值yp（δ），从而可以保证后续数值计算的收敛阶不降低.另外，ur（x）与up（x）相比逼近精度更高，收敛域更大.

取δ=0.4，T=20，在区间[δ，T]上分别使用分段线性（m=1）和分段二次（m=2，d1=0.5，d2=1）插值求数值解.由图2 可知，当x≤2 时ur（x）仍具有O（10-8）的精度，故而在区间[δ，2]上可利用其验证相应配置解un，i（i=1，2，…，m）的正确性，其图像见图3.由图3 可见，配置解和Padé 逼近解在初始阶段的图像完全重合，但随着x 的增加，它们的图像逐渐分开，Padé 逼近解不再有效，但配置解仍然正确，说明方程（39）在大区间上的数值解仍是稳定的.

图3 例2 在大区间上的数值解图像Fig.3 Plots of solutions on large interval in Example 2

5 结语

针对两项线性分数阶微分方程初值问题设计了一种具有较高收敛阶的数值方法.首先将微分方程转化为等价的积分方程，利用Picard 迭代求出解在零点的Puiseux 级数展开式，基于其局部逼近精度高的特点，将原积分方程的奇性分离出来；然后利用得到的级数解的Padé 逼近在小区间上做高精度数值积分计算，保证了整体算法的稳定性，在剩余区间上利用分段Lagrange 插值构造配置方法，并分析了配置格式的收敛性.数值算例结果表明，与在全区间上直接使用配置法相比，所提方法有效提高了数值解的收敛阶.