非線性分數階微分方程的hp 型Legendre 譜配置法

李 珊，安 筱，孫桂磊

（上海理工大學 理學院，上海 200093）

1 問題的提出

分數階微分方程最重要的是非局部性質,它能有效地避免整數階導數的局部性,所以，分數階導數可以更加精確地描述對歷史有依賴的問題。這種非局部性質使得分數階微分方程在許多工程和科學領域，特別是在描述具有記憶過程、遺傳特性[1-3]和異質材料等方面越來越有吸引力。分數階微分方程已被廣泛應用于各種復雜的黏性流體流動模型，如異常擴散過程[4]、分數動力學[5]、粘彈性材料[6]、輸送管道邊界層效應[7]、氧氣通過毛細血管輸送到組織[8]等。同時，利用分數階導數描述問題時，使用少量的參數就能獲得與實際相吻合的結果，所以，與整數階微分方程相比較，在處理復雜系統時，使用分數階微分方程建模將更加簡單、參數意義更加清楚。一般來說，對分數階微分方程的求解是非常困難的。近幾年，學者們對分數階積分微分方程的基本理論及其數值求解進行了大量的研究。因為，一方面分數階微分方程的求解可以轉化為等價的積分方程進行求解；另一方面為使得數學模型更加符合客觀實際，對分數階積分微分方程的求解已經不可避免，所以，對分數階積分微分方程的理論及其數值方法的研究越來越受到學者們的關注。本文提出了一種求解以下非線性分數階初值問題的有效數值方法。

其中，α∈(0,1),f∈C(0,T),并且 Dα是根據下面定義的 α階Caputo 導數。

Caputo 導數 Dα不同于傳統的微分算子，它是一種全局算子。目前分數階導數主要分為3 種類型：Riesz 意義下的分數階導數、Riemann-Liouville分數階導數和Caputo 型分數階導數。1832 年Liouville 采用級數的形式給出了分數階導數，1853年Riemann 采用定積分的形式給出了另一種分數階微分，Grünwald 和Krug 將Riemann 和Liouville的結論進行了統一，得到了Riemann-Liouville 分數階導數的定義。1967 年Caputo 提出了Caputo 型分數階導數。Caputo 型分數階導數因其廣泛的應用背景和相對簡便的運算被廣泛地應用于各個領域。目前在分數階微分方程理論的研究中，關于其數值解進行了大量的研究。在Caputo 型分數階定義下的初值條件具有明確的物理意義，因而在工程物理學領域該方程具有重要的應用價值。但由于解的復雜性，不能直接應用于實際工程中，所以，求解Caputo 型分數階微分方程的數值解受到廣泛關注。

譜方法作為一種全球性的高精度數值方法，在求解微分方程過程中扮演著重要的角色，非常適用于求解分數階微分方程。然而，現有的譜方法主要是一步法。這種方法要求解具有高度正則性。然而，通常對于分數階問題，解不可能是任意光滑的，即使源項非常光滑。因此，一步法缺乏有效處理局部弱奇異的能力。最近，Guo 等[9]提出了一種求解單分數階邊值問題的hp 型Chebyshev譜配置方法，并從理論上和數值上得到了hp 收斂性。在他們工作的基礎上，本文提出了一種求解非線性分數階初值問題的hp 型Legendre 譜配置法，該方法靈活地采用時間步長和局部逼近階數，可以實現任意精度的數值解。并且，本文所提出的數值格式對于數值實現和誤差分析來說更加簡單易行。

2 預備知識

設Ih是區間I=[0,T]上的一個網格,Ih:={tn:0=t0＜t1＜···＜tN=T}，并且設hn=tn-tn-1，In=(tn-1,tn]。通常情況下，對于一個給定區間 Λ和某一權函數 χ(x)，定義

設Lk(x),x∈(-1,1)是標準的階數為k的Legendre多項式。Legendre 多項式集合構成一個完備的L2(-1,1)正交系統，即

其中，δk,j是Kronecker 符號。

定義k階變換的Legendre 多項式

Ln,k(t),k≥0的集合構成一個完備的L2(In)正交系統，即

根據標準Legendre-Gauss 求積的性質，可得

正如文獻[10]中提到的，由于弱奇點的存在(e.g.,ν∈(0,1))，考慮到加權插值求積公式的權值取決于核中的弱奇異因子(t-s)ν是很自然的,對于任意的 φ(s)∈,引入加權插值求積公式

3 hp 型譜配置法

方程(1)可以轉換為以下等價的Volterra 積分方程[12]：

本文的主要目的是基于等價式(13) 提出并分析一種有效的求解方程(1)的 hp型Legendre-Gauss譜配置法。

設un(t)是第n個區間上式(13)的解，即

un(t):=u(t),t∈In,1≤n≤N

由式(13)可得，對于任意t∈In，

求解式(14)的hp 型譜配置法的目的是尋找Un(t)∈，使得

其中，Uk(t)∈是區間Ik上uk(t)的數值解。用U表示全局數值解,即

U(t)=Uk(t),t∈Ik,1≤k≤N

現在描述數值實現,并給出了式(15)的算法。設

由式(6)和式(9)可得，對于 1≤k≤n-1,

其中，

由式(6)，(8)和(11)可得

其中，

由式(15)～(20)可得

根據 Legendre 多項式的正交性，比較式(11)的展開系數得到，對于 0≤p≤Mn，

可采用迭代法計算系統式（22）的展開系數（例如，Newton-Raphson 迭代）。

4 數值結果

現提供一系列數值實驗來驗證所提出算法的有效性。首先，定義在網格點上離散的L2誤差和L∞誤差如下：

例1考慮問題

問題(23)的真解是u(t)=tα。

在圖1（a）和1（b）的第1 個圖中，分別列出了問題（23）在相同的M（多項式階數）以及變化的h（網格長度）下的數值誤差（h型）和，其中，選取M=1，2，4，6。由此可以觀察到h 型的快速收斂速度。

在圖1(a)和1(b)的第2 個圖中，分別列出了問題(23)在相同的h（網格長度）以及變化的M（多項式階數）下的數值誤差，其中，選取h=1，，由此可以觀察到數值誤差的收斂速度隨著M的減小按代數衰減。