重心插值配点法及其应用

重心插值配点法及其应用
摘要：重心Lagrange插值具有数值稳定性好、计算精度高的优点。

采用重心Lagrange插值近似未知函数，建立未知函数各阶导数的微分矩阵。

采用微分矩阵近似未知函数的导数，利用配点法将控制方程和边界条件离散为代数方程组，通过求解代数方程组，从而可求解偏微分方程。

数值算例表明，重心插值配点法具有原理简单，易于程序实现和数值计算精度高的优点。

关键词：重心Lagrange插值；微分矩阵；配点法
Barycentric Lagrange interpolation collocation Method and its Application
Abstract：Barycentric Lagrange interpolation collocation method has excellent numerical stability and high accuracy. this paper presents the Barycentric Lagrange interpolation collocation method to get the differentiation Matrix of unknown function. So the control equation can be expressed as linear systems by the collocation method. According to those formulas, differential equations can be soluted. The principle of this method is simple and easy to programming. The accuracy and the numerical stability are very excellent.
Key words: Barycentric Lagrange interpolation, differentiation Matrix, collocation method,
0 引言
具有初、边值条件的常、偏微分方程的解析解常无法通过理论推导获得，一种有效的途径是采用数值求解方法获得具有一定数值精度的近似解，这些数值方法包括：有限差分法、有限单元法、边界单元法、无网格法及一些杂交使用的方法等数值求解方法。

其中，有限差分法、有限单元法这两种方法要对求解区域划分单元，计算精度依赖于单元的大小。

采用配点法求解边值问题不需要划分单元，公式简单，不需要积分，易于编程。

目前用于求解工程中的常微分方程边值问题的配点法主要有拟谱法和微分求积法。

拟谱法是根据谱方法发展出来的一种方法，虽然这种方法的理论研究已经有进一步的发展，但是工程技术人员对这种方法不是很了解。

微分求积法的基本原理是将未知函数在区间上所有离散点的函数值的加权和来逼近该函数在某一离散点的偏导数或者积分，这种方
法中的权系数的确定通常是根据Lagrange多项式在网点处的导数值给出。

但是这种方法的局限性是离散点不能取得太多，否则Lagrange多项式表示的曲线随多项式次数的升高而出现Runge现象，从而产生计算的不稳定性。

重心Lagrange插值具有极好的数值稳定性和极高的近似精度，同时重心Lagrange插值公式具有紧凑的各阶导数的计算公式。

本文所采用的重心插值配点法就是用重心Lagrange插值多项式求出某一函
数在各个离散点的微分矩阵，从而可以通过矩阵的运算来求出微分方程的解。

而许多结构力学问题的求解最终也归结为在一定的边界条件和初始条件下的(偏)微分方程(组)的求解，所以，可以把重心插值配点法作为用于结构力学问题的求解的一种数值方法。

本文把重心插值配点法用于矩形薄板的挠度问题的分析中，精度较高，得到了满意的效果。

1微分矩阵
设为一待求函数，它在区间[0，L]上连续可微，现沿轴设置个节点，并以节点函数值作为基本未知量，且令=。

在全域内采用重心Lagrange高阶多项式插值逼近，则有：
（1）
其中，为重心权，
（2）
令：（3）
因为所求解的问题是线性的，因此可以定义一个的矩阵。

则（3）式可以写成如下的形式：
（4）
其中为未知函数一阶导数值构成的列向量，为节点函数值的列向量。

就是对求一次导数的微分矩阵。

由（1）、（2）、（3）式可得：
（5）
（6）
同理，可以推得求次导数的微分矩阵。

（7）
（8）
令：（9）
则：（10）
于是，通过矩阵的运算就可以求出节点的函数值，再由（1）式可求出待求
函数各非节点处的值。

2 重心插值配点法求解热传导方程
考虑一维热传导问题，控制微分方程如下：
（11）
其中：为无量纲的时间；为无量纲的轴向坐标；为一个无量纲的参数，现取值为1．上式边界条件为：当=0时，；当=1时，=1．按照重心插值配点法的规则，控制方程及边界条件离散为：
（12）
（13）
（14）
写成矩阵形式如下：
（15）
式中：为的参数列阵；分别是行列的1阶和2阶微分矩阵。

令=,考虑边界条件：则边界约束条件可以表示为：
（16）
将2个方程的左边分别取代矩阵的第1，行，并把矩阵的第行放到矩阵的第1行和第2行之间，则原矩阵的第2行，第3行，第-1行依次成为第3行，第4行，第行。

同样，把第列放到矩阵的第1列和第2列之间，则原矩阵的第2列，第3列，第-1列依次成为第3列，第4列，第列。

将此方程的右边分别取代0矩阵的第1，行。

则根据举证运算即可求得值。

微分矩阵以及节点值的计算，本文都编制Matlab程序来实现。

计算结果及比较如下：
表1温度分布的重心插值解
3、重心插值配点法求解梁弯曲
细长梁在集度为的荷载作用下，在弹性范围内发生线性弯曲，其控制方程为：
（17）
在梁上取个节点，则可以用微分矩阵对个节点值求导表示，写成矩阵形式为：
（18）
式中为的待定的位移参数列阵，是一行列的4阶微分矩阵，为列的载荷列阵。

考虑边界条件：以两端简支梁且在边界有集中力偶为例，当时，。

当时，
则边界约束条件可以表示为：
(19)
将此4个方程的左边分别取代微分矩阵的第1，，2，+1行，将此4个方程的右边分别取代载荷列阵的第1，，2，+1行。

令，则由（18）式可得：
（20）
由（20）式即可求得梁的节点位移列阵。

方程（20）是适定的，可以得到唯一解[9]。

微分矩阵以及梁节点位移的计算，本文都编制Matlab程序来实现。

算例2，两端简支梁，长50cm，高2cm，宽1cm，弹性模量为10Gpa，受如图1所示均布载荷，，得到梁轴线的各节点挠度、转角、弯矩与解析解比较如表1所示（选取5个节点计算）：
图1 受均布载荷梁示意图
表2 梁各点挠度解析解与本文解比较挠度×10-3m
图3 梁挠度示意图
4、结语
重心插值配点法是一种快速求解偏微分方程的方法。

从数值算例中可以看出，这种方法的精度很好，即使是计算复杂分布荷载，也能得到相当高的精度。

在重心插值配点法中采用微分矩阵记号，可以得到简洁的矩阵形式的计算公式，极大地方便了边界条件的离散。

重心插值配点法数学原理简单，编程容易实现，在结构力学的分析中有非常良好的应用前景。

参考文献
王勖成，邵敏. 有限单元法基本原理及数值方法（第2版）[M]. 北京：清华大学出版社，2000.
Weideman, J A C, Reddy, S C. A Matlab differentiation matrix suite[J]. ACM Transactions on thematical Software, 2000, 26(4): 465-519.
Trefethen L N. Spectral methods in Matlab[M]. Philadelphia: SIAM, 2000.
Shu C，Du H. Implementation of clamped and simply supported doundary conditions in the GDQ free vibration analysis of beams and plates[J].Int.J.Solids Structures. 1997,34(7):819-835
聂国隽，仲政. 微分求积单元法在结构工程中的应用[J]. 力学季刊，2005，26（3）：423-427。

Berrut J-P, Trefethen L N. Barycentric Lagrange interpolation[J]. SIAM Review, 2004, 46 (3): 501-517.
王兆清，李淑萍，唐炳涛. 任意连续函数的多项式插值逼近[J]. 山东建筑大学学报,2007,22（2）:158-162.
Schneider C , Werner W. Some new aspects of rational interpolation[J ] . Math Comp , 1986 , 47(175) : 285-299.
李淑萍, 王兆清, 唐炳涛. 重心插值配点法求解初值问题[J]. 山东建筑大学学报, 2007, 22(6):481-485.
赵晓伟，鹿晓阳，汪洪星. 重心插值配点法分析梁弯曲问题[J]. 山东建筑大学学报，2008，23(1)：61-64.
徐芝纶．弹性力学(下册) [M]．北京：高等教育出版社，1990．。