正交配置方法求解一维拟线性抛物方程

王佩臣，张志维，樊玉环，张可为

(黑龙江工程学院 理学院，哈尔滨 150050)

考虑下面一维拟线性抛物方程的初边值问题：

(1)

式中：ψ1(t),ψ2(t),f(u,x,t)和φ(x)为已知函数。

该方程与其他学科有着密切的联系，它的应用直接或间接地源自物理、化学和生物等学科领域，具有非常广泛的实际应用价值。提出一个高阶数值方法求解拟线性抛物方程是一个重要的任务，当前求解该方程主要有有限差分法、有限元、谱方法等，但很少有在时间和空间都比较稳定且有高阶精度的方法。文中提出一个新的方法，组合正交配置方法[1-6]和DIRK[7-8]法来求解一维拟线性抛物方程，该方法在时间和空间都有很高的精度和稳定性。

1 DIRK法

文中简介求解一阶初始值问题的DIRK法，并给出一个具体的四阶A稳定方法。

考虑下面一阶初始值问题：

y′(t)=f(t,y),y(t0)=y0,t≥t0.

(2)

式中：y∈Rm,f:R×Rm→Rm，q步DIRK法可表示为

(3)

tn=t0+nh(n=0,1,…),其中，h>0为积分步长。式(2)、式(3)称为q级龙格库塔法，系数bi,τi,aij(i,j=1,2,…,q)为实数，bi为权系数，τi为节点系数，A=(aij)q×q为系数矩阵。这些系数通常表示为

AτbT

它称为Butcher’s矩阵。取文献[7]中Butcher’s矩阵如下

141416-1121449+416073+124115024-19413001449-4160a41a42a43141015372 091-8794112 1362 091+8794112 136141015372 091-8794112 1362 091+8794112 13614

其中

对于DIRK法，计算ki(i=1,2,3,4)可以用Newton-Raphson法。Newton-Raphson法计算过程如下：

2 正交配置方法

首先，用正交配置方法离散方程(1)的空间导数，得到一个微分方程组，选取Chebyshev-Gauss-Lobatto节点作为配置节点，Legendre多项式作为投影多项式。

设式(1)的解为

(4)

取N个配置点a=x1

(5)

由式(5)求解αi(t)得

(6)

由式(4)得u(x,t)的一阶和二阶偏导为

(7)

(8)

把式(6)代入到式(7)、式(8)中得

(9)

(10)

(11)

(12)

下面选用正交函数作为实验函数,首先考虑下面[-1,1]区间的Legendre多项式。

Legendre多项式：

满足下面性质：

1)递推关系。

2)Legendre多项式的正交性。

其中

下面做y∈[-1,1]→x∈[a,b]的映射。

令[a,b]区间的Chebyshev-Gauss-Lobatto节点

类似上面的讨论，令同样设gi(x)=Pi-1(x)

(13)

式(13)写成多项式展开的形式

(14)

由式(14)得在xj点u(x,t)和u(x,t)的一阶和二阶偏导为

(15)

(16)

(17)

把式(14) 、式(15)、 式(16)写成矩阵的形式为

(18)

其中

U(t)=(u(x1,t),u(x2,t),…,u(xN,t))T,

B(t)=(β1(t),β2(t),…,βN(t))T,

C=(cij)N×N,D=(dij)N×N,

则由式(18)得Β(t)=Q-1U(t)，

(19)

(20)

由式(20) 化简式(1)得

i=2,3,…,N-1.

(21)

其中,u(x0,t)=ψ1(t),u(xN,t)=ψ2(t)，式(21)化为

miNψ2(t)+f(ui,xi,t),i=2,3,…,N-1.

(22)

式(22)满足下面初始条件

u(xk,t0)=φ(xk),k=1,2，…,N-1.

(23)

用第二部分的DIRK法求常微分方程组式(22)、式(23)，进而得到原方程的解。

3 数值实例

3.1 测试实例[9]

初始条件为

u(x,0)=sin (πx), 0

u(0,t)=u(1,t)=0,t≥0.

该问题的精确解为

u(x,t)=e-tsin (πx).

图1为精确解曲面, 图2为在h=0.01,N=11时的数值解曲面，图3为在h=0.01,N=11时的误差曲面。

图1 精确解曲面

图2 在h=0.01,N=11时的数值解曲面

图3 在h=0.01,N=11时的误差曲面

3.2 测试实例

初始条件为

u(x,0)=sinh(πx), 0

u(0,t)=0,u(1,t)=e-tsinh(π),t≥0.

该问题的精确解为

u(x,t)=e-tsinh(πx).

图4为精确解曲面, 图5为在h=0.01,N=11时的数值解曲面，图6为在h=0.01,N=11时的误差曲面。

图4 精确解曲面

图5 在h=0.01,N=11时的数值解曲面

图6 在h=0.01,N=11时的误差曲面

4 结束语

组合使用正交配置法和对角隐式龙格库塔法求解一维热传导方程，使用正交配置法离散空间导数，用对角隐式龙格库塔法求线性解常微分方程组，通过数值实例，可以说明该方法有很高的精度和稳定性，数值实例证实该方法的有效性。