三階Emden-Fowler型微分方程的第4類Chebyshev混合函數數值解法

熊臨晨，許小勇，朱 婷

(東華理工大學理學院，330013，南昌)

0 引言

三階Emden-Fowler型微分方程是與三階常微分方程相關的三階奇異邊值問題，被用來模擬恒星結構建模[1]、人腦中的熱傳導[2]、星系團[3]等研究。Emden-Fowler型微分方程產生于天體物理學中關于氣體動力學的早期研究。從根本上說，它被引入來研究產生于天體物理學中關于氣體動力學的球形氣體云的質量的平衡結構[4]。Emden-Fowler型微分方程由瑞士物理學家J R Emden與統計物理學家R H Fowler共同研究而命名的，具體形式如下[5]：

(1)

受到2類條件約束。

第1類條件：

y(0)=α,y′(0)=β,y″(0)=γ,

(2)

第2類條件：

y(0)=α1,y′(0)=β1,y′(1)=γ1,

(3)

其中α,β,γ,α1,β1,γ1是常數。

近年來，許多學者都在研究三階Emden-Fowler型微分方程。由于這類方程在原點處存在奇異性，求解三階Emden-Fowler型微分方程一直是一個非常具有挑戰性和難度的問題。現有的處理三階奇異問題模型的數值和解析方法很少，如小波方法[6]、樣條方法[7]、變分迭代法[8]、人工神經網絡[9-10]、微分變換法[11]、混合塊法[12]等。

本文提出了一種新的求解三階Emden-Fowler型微分方程的數值方法。該方法是基于第4類Chebyshev多項式和BPF函數的結合，稱為第4類Chebyshev混合函數。由BPF函數與Legendre多項式[13]、Taylor多項式[14]、Lagrange多項式或Bernoulli多項式[15]組合組成的混合函數已被許多學者用來解決數學問題。已有學者提出了第1類[16]、第2類[17]和第3類[18]Chebyshev混合函數，因此本文基于分數階積分公式，引入第4類Chebyshev混合函數，結合配置法將三階Emden-Fowler型微分方程轉化為代數方程組，然后應用牛頓迭代法進行數值求解。

1 基礎知識

分數階導數和積分有很多種定義，其中廣泛使用的分數階導數定義是Caputo導數，分數階積分定義是Rieman-Liouville積分。

定義1：Caputo分數階微分定義為：

(4)

Caputo型分數階微分部分性質：

2)Dα(Iαf(t))=f(t),

3)Dα(λ1f1(t)+λ2f2(t))=λ1Dαf1(t)+λ2Dαf2(t)。

定義2：Rieman-Liouville分數階積分定義為：

(5)

其中*表示tα-1與f(t)的卷積。

定義3：局部可積函數f(t)的Laplace變換F(s)定義如下：

(6)

其中s為復數，該運算具有如下性質：

1)L[λ1f1(t)+λ2f2(t)]=λ1L[f1(t)]+λ2L[f2(t)],

2)L[f(t-α)u(t-α)]=e-αsF(s),

3)L[f*g]=L[f(t)]*L[g(t)]。

定義4：BPF函數是定義在區間[0,T)上的一組階梯函數，滿足

(7)

其中n=1,2,...,m′，對于任意的ψi(x),ψj(x)滿足

2 混合函數及其性質

定義5：第4類Chebyshev多項式[19]Wm(x)是定義在[-1,1]上的n次正交多項式：

(8)

其中x=cosθ，這些正交多項式Wm(x)滿足下列性質

(9)

Wm+1(x)=2xWm(x)-Wm-1(x),m=1,2,...,

初始值為W0(x)=1,W1(x)=2x+1,W2(x)=4x2+2x-1。

第4類Chebyshev多項式的解析形式[20]可以表示為：

(10)

2.1 混合BPF函數與第4類Chebyshev多項式

混合函數φnm(t),n=1,2,...,N,m=0,1,...,M，在區間[0,1)上的定義為

(11)

其中n和m分別是BPF函數和第4類Chebyshev多項式的階數。

2.2 函數逼近

由于混合函數的完備性，任意函數f(t)∈L2[0,1]能被混合函數展開：

(12)

其中

C=[c10,c11,...,c1(M-1),c20,c21,...,c2(M-1),...,cN0,cN1,...,cN(M-1)]T,

Ψ(t)=[φ10,φ11,...,φ1(M-1),φ20,φ21,...,φ2(M-1),...,φN0,φN1,...,φN(M-1)]T。

3 第4類Chebyshev混合函數的Rieman-Liouville型分數階積分公式

現在將積分運算Iα作用于Ψ(t)，有

IαΨ(t)=Ψ*(t)

(13)

其中Ψ*(t)是NM×1維向量，表達式如下

Ψ*(t)=[Iαφ10,Iαφ11,...,Iαφ1(M-1),Iαφ20,Iαφ21,...,Iαφ2(M-1),...,IαφN0,IαφN1,...,IαφN(M-1)]T，

為了得到Iaφnm，使用Laplace變換。

(14)

其中

(15)

對φnm(t)進行Laplace變換

把方程(10)代入上式得到

由Laplace變換性質3)得到

進行Laplace逆變換

(16)

由式(15)和式(16)可得

Iaφnm(t)=

4 誤差與收斂性分析

本小節中證明了第4類Chebyshev混合函數在L2范數意義下的收斂性，為了完成收斂性的證明我們先推導出混合函數展開的誤差上界。

定理1：假設f(t)∈L2[0,1]有一個有界的二階導數，其中|f''(t)|≤B，那么就有以下誤差上界：

(17)

證明：根據L2[0,1]空間的L2范數定義有

(18)

令2Nt-2n+1=x，則

把x=cosθ代入上式，結合Chebyshev多項式定義

記

θsinθdθ。

即上式變為

(19)

接下來分別對I1和I2進行分析計算，首先

由柯西施瓦茲不等式

由sinθ≤1，得

(20)

(21)

由式(19)得

由式(18)得

根據積分判別法，得

則

(23)

定理得證。

5 算法描述

現在給出兩類不同條件下Emden-Fowler型微分方程的混合函數配置點法。

5.1 第1類條件

設方程(1)中未知函數的三階導數可由混合函數表示，即

y?(t)=CTΨ(t),

(24)

對方程(24)從0～t分別進行3次積分并結合條件(2)，可得

y″(t)=CTIΨ(t)+γ,

(25)

y′(t)=CTI2Ψ(t)+tγ+β,

(26)

(27)

將式(25)、式(26)、式(27)代入方程(1)中，得

(28)

5.2 第2類條件

設方程(1)中未知函數的三階導數可由混合函數表示，即

y?(t)=CTΨ(t),

(29)

對方程(29)從0～t分別進行3次積分并結合條件(3)，可得

y″(t)=CTIΨ(t)+y″(0),

(30)

y′(t)=CTI2Ψ(t)+ty″(0)+β1,

(31)

(32)

把y′(1)=γ1代入式(31)中

y″(0)=γ1-CTI2Ψ(1)-β1,

(33)

把式(33)先代入式(30)、式(31)、式(32)中，然后把表達式代入方程(1)中，得

(34)

通過配點法與牛頓迭代法，得到了第4類Chebyshev混合函數的系數cnm，把系數代入式(32)，就能解出方程的數值解。迭代初始值的選取，本文先移除方程中的非線性項得到線性方程，將所求線性方程的解當作迭代初始值。

6 算例分析

例1：考慮帶第1類條件的Emden-Fowler方程

y(0)=0,y′(0)=0,y″(0)=0,

方程的精確解為y(t)=ln(1+t3)。

表1給出了本文方法在不同情況下精確解與數值解的比較，由表1可知絕對誤差隨混合函數參數N、M的提高，越來越小。

表1 不同情況下的精確解與數值解比較

例2：考慮帶第2類條件的Emden-Fowler方程

y?(t)-ty(t)=et(t3-2t2-5t-3), 0

y(0)=0,y′(0)=1,y′(1)=-e。

方程精確解為y(t)=t(1-t)et。

圖1給出了本文的方法在M=8時，N分別為 2、3、4和5的計算結果。從圖1中可以看出，當固定M而提高分辨率尺度N時，方程數值解與精確解之間的誤差隨N增加而變小。表2給出了本文算法與文獻[6]中Bernoulli小波(BWCM)和Hermit小波(HWCM)兩類小波算法在不同參數下絕對誤差的數據比較。不難看出本文方法的誤差更小，精度更高。

圖1 不同參數情況下絕對誤差的數值對比的半對數圖

表2 不同參數情況下本文方法與兩類小波方法的誤差

例3：考慮帶第2類條件的Emden-Fowler方程

y(0)=0,y′(0)=0,y′(1)=4e,

方程精確解為y(t)=t3et。

表3給出本文方法的數值解與精確解的比較。由表3可以看出，在固定分辨率尺度N=8的值，通過提高參數M，精確解與數值解之間絕對誤差越來越小，可以得到更加有效的精確解。

表3 固定N,不同M值情況下的精確解與數值解比較

7 小結與展望

本文提出了一類求Emden-Fowler型微分方程的近似解的Chebyshev混合方法。通過第4類Chebyshev多項式的解析形式與Rieman-Liouville分數積分定義，推導出第4類Chebyshev混合函數的分數階積分公式。利用所推導的積分公式結合配置法將Emden-Fowler型微分方程轉化為一組非線性代數方程，然后應用牛頓迭代法進行求解，降低了方程的復雜程度，并通過數值算例分析證明了算法的有效性與可行性。