移动最小二乘法导数近似讨论

杨建军 杨子乐 黄旺 文丕华

摘要：在移动最小二乘法（moving least squares method， MLS）构造无网格形函数的数值方法中，通常采用无单元伽辽金法（element-free Galerkin method， EFG）的建议，将系数向量a参与导数运算。为探讨这种导数近似算法在更一般无网格法中的适用性和合理性，针对系数向量a是否应参与运算的问题进行讨论和数值检验。结果表明：单纯从近似意义上讲，这种将系数向量代入导数运算的算法并不具有优势；从数值方法的应用意义上讲，这种导数近似算法对数值求解，特别是强式无网格法，会带来一系列潜在不稳定的问题。建议在MLS导数近似中，系数向量a不应当参与导数运算，并提出采用一种由核基函数代替普通基函数的核近似法。

关键词：无网格法； 移动最小二乘法； 导数近似； 系数向量； 核近似

中图分类号：O241.3，O242.2，O302

文獻标志码：A

文章编号：1006-0871（2018）01-0028-07

Abstract： In the numerical method of meshless shape function constructing with moving least squares method（MLS）， the element-free Galerkin method（EFG） is adopted customarily， and the coefficient vector a usually take part in the derivative operation. To investigate the applicability and rationality of this derivative approximation algorithm applying to more general meshless method， the discussion and numerical tests are carried out which focus on whether the coefficient vector a should take part in computing or not. The results show that on approximate means only， there is no advantage if the coefficient vector takes part in computing； on the application means of the numerical method， there is a series of potential instability problems using the derivative approximation algorithm on numerical calculation， especially for the strong meshless method. It is suggested that the coefficient vector a should not be included in the derivative operation. A core approximation method is proposed in which the common basis function is replaced by core primary function.

Key words： meshless method； moving least squares method； derivative approximation； coefficient vector； core approximation

0 引 言

移动最小二乘 （moving least squares， MLS）法最初用于地理信息绘制和表面拟合。[1-3]随着无网格法研究的兴起，MLS法被广泛用于构造无网格形函数。[4-9]

BELYTSCHKO等[10]提出无单元伽辽金法（element-free Galerkin methods， EFG）时，主张系数向量a（x）不应当被假定为常数向量，而应参与求导运算。此后，这一观点被广为接受，在MLS法的导数近似中得到普遍应用。[11-18]

MLS法在执行导数近似时，系数向量a（x）是否该参与导数运算是一个由来已久的问题。鉴于该问题涉及导数近似中基本概念的争议，并且应用范围非常广泛，因此有必要重新审视。

1 移动最小二乘法

设问题域Ω上的场函数为u（x），这个场函数的值只在一些域内的散乱点上已知或可测量。为获得目标点（计算点）

xI处的场函数uI（x），为xI设置一个以其为中心的局部支撑域ΩS，见图1。只要这个支撑域足够局部，场函数u（x）总可以用一个预定义的试函数

u～（x）近似描述。

由于多项式函数较简单，因此MLS法一般采用多项式试函数，即

以上所述即为经典的MLS法。在此基础上，可以使用复变量近似方法提高计算效率和计算精度[21]，或者通过改造基函数使其具有插值性[22-23]。此外，为提高MLS法的计算稳定性和计算精度，可使用核基函数代替普通基函数，即移动最小二乘核近似（moving least squares core，MLSc）方法。[20，24-25]式（2）对应的核基函数可写为

MLSc法可以保证A矩阵条件良好，从而可以显著提高近似计算的精度和稳定性。[20]这一优势将在后续数值算例中给出检验。

2 MLS导数近似

式（12）是目标点xI处场函数的近似，在数值方法中，如果要执行场函数导数的近似，可通过形函数J的求导实现，即

上述2种导数近似方法的本质区别是向量a作常数考虑还是作非常数考虑，对该分歧讨论如下。

（1）式（4）引入权函数w（x），但不应当改变向量a的初始假定和本质属性，a作常数考虑合理。

（2）a不作常数考虑，则除基函数需求导外，还要对权函数求导，无疑会增加对权函数连续性的要求。通常，权函数的连续性要高于基函数的连续性，权函数求导并不会提高近似函数的连续性，因为近似函数的连续性由基函数的连续性决定。

（3）a不作常数考虑时，形函数J的导数中会出现A及其导数求逆的累加，如式（21）和（22）。从原理上讲，这会造成误差累计，对计算精度不利。

（4）MLS法本质上是一种非常简单的方法，但a参与导数运算后，使算法变得非常复杂，這样会削弱MLS法的竞争力。

（5）如果a不作常数考虑确定对EFG法[10]有利，那么对更一般的方法，这种优势是存在疑问的。

接下来，进一步通过数值试验来进行审视。

3 数值检验

3.1 算例1

曲面拟合问题对近似方法的检验比较直接，因此通常被用于近似算法检验。假设在x∈[0，1]且y∈[0，1]问题域上的曲面函数为

其函数图像见图2。该函数的偏导数很容易得到，此处从略。在近似计算中，为避免散乱节点的扰动，采用规则节点离散方案。节点间距h=0.2的场节点离散方案见图3。

在近似计算中，任一目标节点xI的场函数u，一阶导数u，x和二阶导数u，xx由其他节点的场函数值通过近似计算得到。

为评价近似计算的精度和收敛性，定义平均L2误差范数为

式中：ξI，num为数值计算得到的值；ξI，exa为解析得到的精确值；N为考察的计算点数量。对于曲面拟合问题，N为所有的场节点数量；对无网格法求解问题，N为结果输出路径上的取值点数量。

设置一组疏密不同的节点离散方案（节点间距为h），将原函数、一阶导数和二阶导数的近似误差汇总，见图4。图例中，“MLS”表示使用经典MLS法，“MLSc” 表示使用改进后的MLSc法，a表示将a视为常数不参与导数运算，a（x）表示将a视为非常数参与导数运算。

忽略一些细节问题，图4可以解读为2个主要结论：（1）a参与导数运算并没有表现出预期的优势；（2）MLS法在节点稠密时收敛性将变差，而MLSc法总是表现为严格的收敛性。

3.2 算例2

为不失一般性，考虑一个形状较为复杂的，或者场函数光滑性较弱的曲面拟合问题。假设在x∈[-8，8]且y∈[-8，8]问题域上的曲面函数为

该曲面函数的图像见图5。

绘制一组疏密不同的节点离散方案和场函数及其导数的近似计算误差，见图6。由此得出的结论与算例1保持一致：（1）a参与导数运算在近似精度和收敛性上，比a不参与导数运算更差；（2）MLS法随节点加密会丧失收敛性，MLSc法总是表现为更高的计算精度和严格的收敛性。

3.3 算例3

结合一种具体的无网格方法进一步检验，使用有限点法（finite point method，FPM）[26-27]。选择FPM主要考虑到：该法是直接的配点法，属于最简单的无网格法之一，对近似计算的验证受数值离散算法的影响较小[28]；更重要是，FPM属于强式法，会使用到二阶导数近似，对问题的验证会更完整。

选择一个常用的悬臂梁问题算例，其结构模型见图7。计算采用无量纲化（量纲归一化）处理，计算参数取：梁宽D=2，梁长L=12，梁端荷载集度P=6，弹性模量E=10-4，泊松比ν=1/3。将y=0路径上的纵向位移uy和x=L/2路径上的剪切应力σxy作为计算指标，其解析解[10，29]为

式中：I为截面惯性矩，I=D-3/12。数值计算中为避免散乱节点扰动，采用规则节点离散方案，但各坐标方向的标准节点间距不同，且hx=2hy。hy=0.4的场节点离散方案见图8。为便于读者直观了解FPM求解该问题的效果，给出hy=0.2离散方案的计算结果，见图9。

2种导数近似方法分别采用MLS法和MLSc法的数值计算收敛性比较见图10（最密节点间距取hy=0.05）。由此可以看出，a不作常数考虑时，不论是位移解还是应力解，其计算收敛性要比另外的方案明显变差，只是在节点较为稀疏时，似乎a不作常数考虑的解更精确。但须注意到，在节点非常稀疏的情况下，求解精度本身是非常粗糙的，不具有代表性，所以此处的结论依然应该是：a不作常数考虑时，求解收敛性并不具有优势。

无网格方法对局部支撑域尺度参数αS的敏感性，不仅直接反映数值方法计算的稳定性，也可以反映近似计算的稳定性。因此，对支撑域尺度参数的数值计算影响效果进行分析，hy=0.2时的计算结果见图11。显然，不论是位移解还是应力解，a参与导数计算的数值结果明显要差得多，换句话说，a参与导数计算的情况下，数值方法的计算稳定性会变差。必须注意的是，支撑域尺度参数αS=3.4时，a参与导数计算的效果出现偶发性的好转，而图9的结果正好采用这一参数，均不具有普遍性，因此此处的结论依然类似，a不作为常数考虑或a参与导数计算时，在数值计算的稳定性上不具有优势。

3.4 算例4

在矩形域Ω（0≤x≤a，-b/2≤y≤b/2）上，当取a=b=1时，泊松方程可定义为

其函数图像见图12。采用FPM求解该问题，问题域采用不规则节点离散，见图13。

取x=0.8的输出路径，选择21个等间距分布的取值点，此结果输出路径上的数值计算收敛性见图14。由此可以得到与前几个算例类似的结论：（1）总体上，a参与导数计算的求解结果会更差，虽然在随机节点高度稠密区间，即ha<0.02时，a参与导数计算的结果似乎更好些，但这个区间对经典MLS而言，数值解已经丧失精确性，其优劣性的评价意义不大；（2）改进的MLSc法的数值解严格收敛，而MLS会随着节点稠密变得非常不稳定。

4 结 论

通过4个数值算例验证，均得出a参与导数运算的负面影响，这显然与EFG方法的建议不同。其中，2个算例是通过曲面拟合的方式直接检验MLS的导数近似效果，对问题的评价具有普遍意义。此外，采用FPM通过2个数值求解算例评价a参与导数运算的优劣性。单用一种具体的方法评价，虽然有失全面性和公允性，但至少可以说明一个问题：EFG方法中认为a参与导数运算使数值求解更精确，这个结论不是对任何方法都适用的，不具有普遍性。

在EFG方法[10]中，对a参与导数运算的效果评价通过一个简单小片试验算例进行验证，场节点非常稀疏，精确解又是一个简单的二次多项式，加上EFG方法只会使用一阶导数近似，因此其论据支持显然偏弱，说服力不够，将其结论不加鉴别的普遍推广，是不恰当的。

从原理性分析的角度来看，系数向量a不参与导数运算是完全合理的。简单地说，多项式函数的求导与系数无关，将a作为函数处理并参与导数运算，会增加问题的复杂性。更重要的是，经过本文的初步数值验证，这种处理会使情况变差而非向好。

此外，数值算例也对核近似方法MLSc和经典的MLS法进行比较，进一步证明MLSc法是一种严格收敛和计算稳定的有效方法。此数值检验结果完全符合文献[20]的理论解释。MLSc法的改进措施很简单，但对提升MLS方法的计算效果却非常显著。

在MLS导数近似中，系数向量a不应当参与导数运算，建议采用更简单的MLSc方法使计算更加精确和稳定，其在无网格近似方法的选择比较中更具有竞争力。

参考文献：

[1] MCLAIN D H. Drawing contours from arbitrary data points[J]. Computer Journal， 1974， 17（4）： 318-324.

[2] GORDON W J， WIXOM J A. Shepards method of “metric interpolation” to bivariate and multivariate interpolation[J]. Mathematics of Computation， 1978， 32（141）： 253-264.

[3] LANCASTER P， SALKAUSKAS K. Surfaces generated by moving least squares methods[J]. Mathematics of Computation， 1981， 37（155）： 141-158.

[4] BELYTSCHKO T， KRONGAUZ Y， ORGAN D， et al. Meshless method： An overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering， 1996， 139（1-4）： 3-47.

[5] 張雄， 刘岩， 马上. 无网格法的理论及应用[J]. 力学进展， 2009， 39（1）： 1-36.

[6] 程玉民. 移动最小二乘法研究进展与述评[J]. 计算机辅助工程， 2009， 18（2）： 5-11. DOI： 10.3969/j.issn.1006-0871.2009.02.003.

[7] YANG J J， ZHENG J L. Intervention-point principle of meshless method[J]. Chinese Science Bulletin， 2013， 58（4-5）： 478-485.

[8] 顾元通， 丁桦. 无网格法及其最新进展[J]. 力学进展， 2005， 35（3）： 323-337.

[9] NAYROLES B， TOUZOT G， VILLON P. Generalizing finite element method： Diffuse approximation and diffuse elements[J]. Computational Mechanics， 1992， 10（5）： 307-318.

[10] BELYTSCHKO T， LU Y Y， GU L. Element-free Galerkin methods[J]. International Journal for Numerical Methods in Engineering， 1994， 37（2）： 229-256.

[11] ATLURI S N， ZHU T. A new meshless local Petrov-Galerkin（MLPG） approach in computational mechanics[J]. Computational Mechanics， 1998， 22（2）： 117-127.

[12] LIU G R. Meshfree methods： Moving beyond the finite element method[M]. 2nd ed. Boca Raton： CRC Press， 2009： 237-272.

[13] ZHU T， ZHANG J D， ATLURI S N. A local boundary integral equation（LBIE） method in computational mechanics， and a meshless discretization approach[J]. Computational Mechanics， 1998， 21（3）： 223-235.

[14] 张雄， 胡炜， 潘小飞， 等. 加权最小二乘无网格法[J]. 力学学报， 2003， 35（4）： 425-431.

[15] WEN P H， ALIABADI M H. An improved meshless collocation method for elastostatic and elastodynamic problems[J]. International Journal for Numerical Methods in Biomedical Engineering， 2008， 24（8）： 635-651. DOI： 10.1002/CNM.977.

[16] 张雄， 刘岩. 无网格法[M]. 北京： 清华大学出版社， 2004.

[17] 陈文， 傅卓佳， 魏星. 科学与工程计算中的径向基函数方法[M]. 北京： 科学出版社， 2014.

[18] 程玉民. 无网格方法[M]. 北京： 科学出版社， 2015.

[19] 左传伟， 聂玉峰， 赵美玲. 移动最小二乘方法中影响半径的选取[J]. 工程数学学报， 2005， 22（5）： 833-838.

[20] 杨建军， 郑健龙. 移动最小二乘法的近似稳定性[J]. 应用数学学报， 2012， 35（4）： 637-648.

[21] 程玉民， 彭妙娟， 李九紅. 复变量移动最小二乘法及其应用[J]. 力学学报， 2005， 37（6）： 719-723.

[22] REN H P， CHENG Y M， ZHANG W. An interpolating boundary element-free method（IBEFM） for elasticity problems[J]. Science China： Physics， Mechanics and Astronomy， 2010， 53（4）： 758-766.

[23] WANG J F， SUN F X， CHENG Y M， et al. Error estimates for the interpolating moving least-squares method[J]. Applied Mathematics and Computation， 2014， 245： 321-342.

[24] 杨建军， 郑健龙. 无网格介点法： 一种具有h-p-d适应性的无网格法[J]. 应用数学和力学， 2016， 37（10）： 1013-1025.

[25] 杨建军， 郑健龙. 无网格局部强弱法求解不规则域问题[J]. 力学学报， 2017， 49（3）： 659-666. DOI： 10.6052/0459-1879-16-383.

[26] OATE E， IDELSOHN S， ZIENKIEWICZ O C， et al. A finite point method in computational mechanics： Applications to convective transport and fluid flow[J]. International Journal for Numerical Methods in Engineering， 1996， 39（22）： 3839-3866. DOI： 10.1002/（SICI）1097-0207（19961130）39：22<3839：：AID-NME27>3.0.CO；2-R.

[27] OATE E， PERAZZO F， MIQUEL J. A finite point method for elasticity problems[J].Computers & Structures， 2001， 79（22-25）： 2151-2163.

[28] 王冰冰， 段庆林， 李锡夔， 等. Euler梁弯曲分析的无网格高阶曲率光顺方案[J]. 计算机辅助工程， 2017， 26（4）： 1-6. DOI： 10.13340/j.cae.2017.04.001.

[29] TIMOSHENKO S P， GOODIER J N. Theory of elasticity[M]. 3rd ed. New York： McGraw-Hill， 2004.

[30] 杨建军， 郑健龙. 无网格全局介点法[J]. 应用力学学报， 2017， 34（5）： 956-962.

（编辑 武晓英）