熊臨晨，許小勇，朱 婷

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

0 引言

三階Emden-Fowler型微分方程是與三階常微分方程相關(guān)的三階奇異邊值問(wèn)題，被用來(lái)模擬恒星結(jié)構(gòu)建模[1]、人腦中的熱傳導(dǎo)[2]、星系團(tuán)[3]等研究。Emden-Fowler型微分方程產(chǎn)生于天體物理學(xué)中關(guān)于氣體動(dòng)力學(xué)的早期研究。從根本上說(shuō)，它被引入來(lái)研究產(chǎn)生于天體物理學(xué)中關(guān)于氣體動(dòng)力學(xué)的球形氣體云的質(zhì)量的平衡結(jié)構(gòu)[4]。Emden-Fowler型微分方程由瑞士物理學(xué)家J R Emden與統(tǒng)計(jì)物理學(xué)家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是常數(shù)。

近年來(lái)，許多學(xué)者都在研究三階Emden-Fowler型微分方程。由于這類方程在原點(diǎn)處存在奇異性，求解三階Emden-Fowler型微分方程一直是一個(gè)非常具有挑戰(zhàn)性和難度的問(wèn)題。現(xiàn)有的處理三階奇異問(wèn)題模型的數(shù)值和解析方法很少，如小波方法[6]、樣條方法[7]、變分迭代法[8]、人工神經(jīng)網(wǎng)絡(luò)[9-10]、微分變換法[11]、混合塊法[12]等。

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

1 基礎(chǔ)知識(shí)

分?jǐn)?shù)階導(dǎo)數(shù)和積分有很多種定義，其中廣泛使用的分?jǐn)?shù)階導(dǎo)數(shù)定義是Caputo導(dǎo)數(shù)，分?jǐn)?shù)階積分定義是Rieman-Liouville積分。

定義1：Caputo分?jǐn)?shù)階微分定義為：

(4)

Caputo型分?jǐn)?shù)階微分部分性質(zhì)：

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

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

定義2：Rieman-Liouville分?jǐn)?shù)階積分定義為：

(5)

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

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

(6)

其中s為復(fù)數(shù)，該運(yùn)算具有如下性質(zhì)：

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函數(shù)是定義在區(qū)間[0,T)上的一組階梯函數(shù)，滿足

(7)

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

2 混合函數(shù)及其性質(zhì)

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

(8)

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

(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多項(xiàng)式的解析形式[20]可以表示為：

(10)

2.1 混合BPF函數(shù)與第4類Chebyshev多項(xiàng)式

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

(11)

其中n和m分別是BPF函數(shù)和第4類Chebyshev多項(xiàng)式的階數(shù)。

2.2 函數(shù)逼近

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

(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混合函數(shù)的Rieman-Liouville型分?jǐn)?shù)階積分公式

現(xiàn)在將積分運(yùn)算Iα作用于Ψ(t)，有

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

(13)

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

Ψ*(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)

對(duì)φnm(t)進(jìn)行Laplace變換

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

由Laplace變換性質(zhì)3)得到

進(jìn)行Laplace逆變換

(16)

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

Iaφnm(t)=

4 誤差與收斂性分析

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

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

(17)

證明：根據(jù)L2[0,1]空間的L2范數(shù)定義有

(18)

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

把x=cosθ代入上式，結(jié)合Chebyshev多項(xiàng)式定義

記

θsinθdθ。

即上式變?yōu)?/p>

(19)

接下來(lái)分別對(duì)I1和I2進(jìn)行分析計(jì)算，首先

由柯西施瓦茲不等式

由sinθ≤1，得

(20)

(21)

由式(19)得

由式(18)得

根據(jù)積分判別法，得

則

(23)

定理得證。

5 算法描述

現(xiàn)在給出兩類不同條件下Emden-Fowler型微分方程的混合函數(shù)配置點(diǎn)法。

5.1 第1類條件

設(shè)方程(1)中未知函數(shù)的三階導(dǎo)數(shù)可由混合函數(shù)表示，即

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

(24)

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

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

(25)

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

(26)

(27)

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

(28)

5.2 第2類條件

設(shè)方程(1)中未知函數(shù)的三階導(dǎo)數(shù)可由混合函數(shù)表示，即

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

(29)

對(duì)方程(29)從0～t分別進(jìn)行3次積分并結(jié)合條件(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)中，然后把表達(dá)式代入方程(1)中，得

(34)

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

6 算例分析

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

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

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

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

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

例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時(shí)，N分別為 2、3、4和5的計(jì)算結(jié)果。從圖1中可以看出，當(dāng)固定M而提高分辨率尺度N時(shí)，方程數(shù)值解與精確解之間的誤差隨N增加而變小。表2給出了本文算法與文獻(xiàn)[6]中Bernoulli小波(BWCM)和Hermit小波(HWCM)兩類小波算法在不同參數(shù)下絕對(duì)誤差的數(shù)據(jù)比較。不難看出本文方法的誤差更小，精度更高。

圖1 不同參數(shù)情況下絕對(duì)誤差的數(shù)值對(duì)比的半對(duì)數(shù)圖

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

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

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

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

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

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

7 小結(jié)與展望

本文提出了一類求Emden-Fowler型微分方程的近似解的Chebyshev混合方法。通過(guò)第4類Chebyshev多項(xiàng)式的解析形式與Rieman-Liouville分?jǐn)?shù)積分定義，推導(dǎo)出第4類Chebyshev混合函數(shù)的分?jǐn)?shù)階積分公式。利用所推導(dǎo)的積分公式結(jié)合配置法將Emden-Fowler型微分方程轉(zhuǎn)化為一組非線性代數(shù)方程，然后應(yīng)用牛頓迭代法進(jìn)行求解，降低了方程的復(fù)雜程度，并通過(guò)數(shù)值算例分析證明了算法的有效性與可行性。