兩項線性分數階微分方程的奇點分離分段配置法

李玉玉，廉 歡，王同科

（天津師范大學數學科學學院，天津 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 插值構造配置方法，并分析了配置格式的收斂性.數值算例結果表明，與在全區間上直接使用配置法相比，所提方法有效提高了數值解的收斂階.