张瑞平， 李成澄

(西安理工大学 理学院， 陕西 西安 710054)



基于7次样条函数的微分求积法研究

张瑞平， 李成澄

(西安理工大学 理学院， 陕西 西安 710054)

针对传统微分求积法的局限性，提出了一种新的微分求积法.该方法采用7次B样条函数构造了区间上的插值基函数，并结合微分求积法的基本思想，形成了基于7次样条的微分求积法.该方法不仅保留了传统微分求积法计算复杂度小、精度高等优点，还具有数值稳定性好，节点可以取很多的优点.数值算例表明了该方法的有效性.

7次B样条函数； 插值基函数； 微分求积法； 对流占优方程

0 引言

微分求积法(differential quadrature method，DQM)最初是由Bellman和casti[1]在1971年根据微积分思想提出的一种求解微分方程的数值方法.1987年，Bert等[2]首次将这种方法应用于结构力学分析，发现该法具有公式简单、使用方便、计算量少、精度高等优点，从而掀起了对DQ法及其在工程中应用的研究热潮.到现在为止，该方法被成功应用到众多领域，如生物学和医学问题、系统识别和系统响应问题、非线性扩散问题、非线性热传导问题、化学反应问题等.传统DQ法有节点不能取多等缺点，于是人们提出了各种改进的DQ法，如调和微分求积法[3]小波—DQ法[4]、移动最小二乘微分求积法[5]、重心插值微分求积法[6]、局部径向基微分求积法[7]等.然而，自然界中许多物理现象都可以用对流占优方程描述，因此构造快速、稳定的数值算法，有着重要的理论和实际意义.针对对流占优问题，许多研究者构造了大量高精度的计算格式来克服数值振荡现象.如特征有限差分格式[8,9]、特征有限元法[9,10]、Streamline upwind Petrov-Galerkin( SU PG) 方法[11]、无网格法[12]等.样条微分求积法最早是由Kashef和Bellman[13]提出的，钟宏志给出了详细的构造方法，形成了5次样条微分求积法[14]，并将其用到求解基尔霍夫板的问题.本文在此基础上构造了7次样条插值基函数，形成7次样条微分求积法，是对B样条微分求积法的有益补充.数值算例表明该方法具有高精度、数值稳定性好的优点.

1 算法构造

首先将正则化区间[0,1]N等分，步长为h，节点为

xj=jh,j=0,1,2,…,N

则7次B样条函数有如下正则化表达式[15]：

(1)

为了在区间上建立插值基函数，需要引入区间以外的点来满足边界条件的要求.这里，我们增加如下节点

x-j=-jh,xN+j=(N+j)h,j=1,2,3,4

建立如下插值函数

(2)

其中pj(x)是要构造的插值基函数,为了能够满足如下插值条件的要求

f7(xi)=yi

插值基函数pj(x)在每个点应该满足如下的基本条件：

(3)

i,j=-4,-3,-2,-1,0,1,…,N,N+1,N+2,N+3,N+4

我们设想用7次B样条函数的组合形式来构造该插值基函数.为了能够达到目的，我们先建立如下的辅助样条函数

(4)

其中

i=0,1,2,3,4

an+1,3(yi-2+yi+2)+an+1,4(yi-3+yi+3)+

an+1,5(yi-4+yi+4),

n=0,1,2,3,4

(5)

其中系数aij的取值如表1所示.

表1 aij(i,j=1,2,3,4,5)的值

由(5)式消去

yi-4,yi-3,yi-2,yi-1,yi+1,yi+2,yi+3,yi+4,

得到如下表达式

(6)

其中qi,i=0,1,2,3,4的值如表2所示.

表2 qi(i=0,1,2,3,4)的值

根据(2)式和(6)式，我们有理由构造如下函数

(7)

再将(4)式代入(7)式，我们发现基函数的表达式如下

j=-4,-3,…,N+3,N+4

(8)

为了消除区间之外的8个节点，我们在区间内两端附近加入如下的8个非整数节点:

x1/5=h/5,x2/5=2h/5,x3/5=3h/5,

x4/5=4h/5,xN-4/5=(N-4/5)h,

xN-3/5=(N-3/5)h,xN-2/5=(N-2/5)h,

xN-1/5=(B-1/5)h,

(9)

将这些非整数节点代入 (2) 式，由插值函数的局部性可得关于y-4,y-3,y-2,y-1的方程组

(10)

和关于yN+1,yN+2,yN+3,yN+4的方程组

(11)

求解方程组(10)和(11)，求得y-4,y-3,y-2,y-1和yN+1,yN+2,yN+3,yN+4的表达式，再代入(2)式重新整理得

(12)

其中φj(x)就是所求的基于7次样条函数u7(x)的插值基函数，具体表达式为

φj(x)=pj(x) 6≤j≤N-6

(13)

这里cij=bi,9-j,i,j=0,1,2,…,9，其中c04=1;cii=1,(i=5,6,7,8,9)，其余cij(i,j=0,1,2,…,9)的值如表3所示.

表3 cij(i,j=0,1,2,…,9)的值

对φi(x)求各阶导数即可得到7次样条微分求积法的各阶权系数.

2 数值算例及分析

以对流占优方程的定解问题为例

其中ε为扩散系数，该问题的精确解为

这里应用基于本文构造的插值基函数的DQ法求解，边界条件用直接法处理，则对流占优方程定解问题可以离散为

(14)

求解方程组即可得到对流占优方程定解问题的数值解，当ε分别取不同值时计算结果如图1所示.

从图1可以看出，应用本文算法的近似解逼近效果很好，在当ε取很小时如(ε=0.005)，该算法表现出了很好的稳定性.

表4给出了本文算法和文献[16]算法计算结果的比较.可以看出，当ε取值较小时，本文算法收敛速度比文献[16]要小，但在取值较大时，本文方法的收敛性优于文献[16].这说明本文方法在解决较平滑性问题时有明显的优势.

(a)ε=1.0对流占优问题的解

(b)ε=0.1对流占优问题的解

(c)ε=0.01对流占优问题的解

(d)ε=0.005对流占优问题的解

表4N取不同值时，不同的ε取值下全部节点的均方根误差

误差ε=1.0ε=0.1ε=0.01ε=0.005N=100本文算法文献[16]算法2.0366×10-121.9312×10-104.2050×10-121.2707×10-92.0381×10-51.0075×10-61.2370×10-31.0303×10-3N=150本文算法文献[16]算法2.0100×10-127.9382×10-112.6688×10-122.3909×10-101.1221×10-62.3084×10-91.2592×10-43.0908×10-5

3 结论

本文以7次B样条函数构造插值基函数，形成了基于7次样条的微分求积法，该方法不仅克服了传统微分求积法的缺点，还继承了样条函数的优点，如稳定性强，灵活性好，精度高等特点，能够有效求解高阶微分方程定解问题.除此之外，7次样条插值基函数形成的权系数矩阵呈稀疏带状，因此为大型工程问题提供了可能.

[1] Bellman R E,Casti J.Differential quadrature and long-term integration[J].Journal of Mathematical Analysis and Applications,1971,4(34):235-238.

[2] Bert C W,Jang S K,Striz A G.Two new approximate methods for analyzing free vibration of structural components[J].American Institute of Aeronautics and Astronautics,1988,26(5):612-618.

[3] Striz A G,Wang X,Bert C W.Harmonic differential quadrature method and appli-cations to structural components[J].Acta Mechanical,1995,111(2):85-94.

[4] 赵凤群，张培茹，张瑞平.两点边值问题的小波配点法[J].计算力学学报，2009,26(6):947-951.

[5] Liew K M,Hang Y Q,Reddy J N.Moving least squares differential quadralture method and its application to the analysis of shear deformable plates[J].International Journal for Numerical Methods in Engineering,2003,56:2 331-2 351.

[6] 王光法.重心插值微分求积法及其在结构工程中的应用[D].济南：山东建筑大学,2009.

[7] 魏学润.弹塑性扭转问题的局部微分求积法[D].苏州：苏州大学,2010.

[8] 秦新强，马逸尘.对流占优扩散方程的一种特征差分算法[J].高等学校计算数学学报,2004,26(1):30-37.

[9] Douglas J J,Russell T F.Numerical method for convection dominared diffusion problem based on combining the method of characteristics with finite element or finite differ ence procedures[J].SIAM Journal on Numerical Analysis,1982,19(5):871-885.

[10] Qin X,Ma Y,Zhang Y.Two-grid method for characteristics finite-element solution of 2D nonlinear convection- dominated diffusion problem[J].Applied Mathematics and Mechanics,2005,26(11):1 365-1 372.

[11] Donea J,Huerta A.Finite element methods for flow problems[M].Chichester:Wiley,2003.

[12] 苏李君，秦新强，王志刚.对流占优扩散方程的无网格法[J].武汉理工大学学报，2010,32(11):172-176.

[13] Bellman R E,Kashef B.Solution of the partial differential equation of the hodgkins-Huxley model using differential quadrature[J].Math.Biosic.,1974,19:1-8.

[14] Zhong H Z.Spline-based differential quadrature for fourth order equations and itsapplication to Kirchhoff plates[J].Applied Mathematical Modeling,2004,28:353-366.

[15] 张锡治，崔旭民.基于七次B样条的结构分析方法[J].工程力学，2009,26(4)：16-21.

[16] 陆小平.微分方程微分求积法研究[D].西安：西安理工大学，2013.

The differential quadrature method based onseven times spline function

ZHANG Rui-ping， LI Cheng-cheng

(School of Science, Xi′an University of Technology, Xi′an 710054, China)

Considering the limitation of the traditional differential quadrature method,a novel different-tial quadrature method is put forward.The interpolation basis function on interval are build by suing seven times B spline function.Cmbined with the basic idea of DQ method, the differential quadrature method based on seven times spline function is put forward.This method not only retains the traditional differential quadrature method which can achieve very high precision, but also has good numerical stability,and the node can take many advantages.Numerical examples show the effectiveness of the proposed method.

seven times B spline function; interpolation basis function; differential quadrature method; convection-dominate equation

2015-01-07

陕西省教育厅专项科研计划项目(12JK0534)

张瑞平(1963-)，男，陕西渭南人，副教授，研究方向:微分方程数值解

1000-5811(2015)03-0177-04

O24

A