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

李 珊，安 筱，孫桂磊

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

1 問題的提出

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

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

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

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

2 預備知識

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

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

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

定義k階變換的Legendre 多項式

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

根據(jù)標準Legendre-Gauss 求積的性質(zhì)，可得

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

3 hp 型譜配置法

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

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

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

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

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

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

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

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

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

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

其中，

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

其中，

由式(15)～(20)可得

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

可采用迭代法計算系統(tǒng)式（22）的展開系數(shù)（例如，Newton-Raphson 迭代）。

4 數(shù)值結(jié)果

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

例1考慮問題

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

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

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