三阶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型微分方程转化为一组非线性代数方程，然后应用牛顿迭代法进行求解，降低了方程的复杂程度，并通过数值算例分析证明了算法的有效性与可行性。