维普资讯 http://www.cqvip.com 第３０卷第２期 南京师大学报（自然科学版） Ｖｏｌ＿３０ Ｎｏ．２ ２００７年６月 ＪＯＵＲＮＡＬ ＯＦ ＮＡＮＪＩＮＧ ＮＯＲＭＡＬ ＵＮＩＶＥＲＳＩＴＹ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ） Ｊｕｎ，２００７ 一种基于变分的网格运动方法及其在 边界层问题数值求解中的应用 李征，王双虎 （北京应用物理与计算数学研究所，北京１０００８８） ［摘要］ 从控制网格以适应解的性态并具有一定光滑性的朴素思想出发，建立了一种基于变分的网格运动方 法，并给出了两种简单的数值求解方法．以两点边值问题及边界层问题为例，对这一网格运动方法进行了数值 试验并得到了满意的数值结果，结果说明本方法成功地控制了网格的分布． ［关键词］ 运动网格，变分法，控制函数，欧拉一拉格朗吕方程 ［中图分类号］０２４２．２１ ［文献标识码］Ａ ［文章编号］１００１－４６１６（２００７）０２－００１５－０７ Ａ Ｍｏｖｉｎｇ Ｍｅｓｈ Ｍｅｔｈｏｄ Ｂａｓｅｄ ｏｎ Ｖａｒｉａｔｉｏｎａｌ Ｐｒｏｂｌｅｍ ｆｏｒ Ｂｏｕｎｄａｒｙ Ｌａｙｅｒ Ｐｒｏｂｌｅｍ Ｌｉ Ｚｈｅｎｇ，Ｗａｎｇ Ｓｈｕａｎｇｈｕ （Ｉｎｓｔｉｔｕｔｅ ｏｆ Ａｐｐｌｉｅｄ Ｐｈｙｓｉｃｓ ａｎｄ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｍａｔｈｅｍａｔｉｃｓ，Ｂｅｉｊｉｎｇ １０００８８，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ａ ｎｅｗ ｍｏｖｉｎｇ ｍｅｓｈ ｍｅｔｈｏｄ ｉｓ ｇｉｖｅｎ ｉｎ ｔｈｉｓ ｄｉｓｓｅｒｔａｔｉｏｎ．Ｋｅｙ ｉｄｅａ ｏｆ ｏｕｒ ｍｅｔｈｏｄ ｉｓ ｔｈａｔ ｔｈｅ ｍｅｓｈ ｓｈｏｕｌｄ ａｄａｐｔ ｔｈｅ ｓｏｌｕｔｉｏｎｓ ｃｈａｒａｃｔｅｒ ａｎｄ ｈａｖｅ ｓｏｍｅ ｓｍｏｏｔｈｎｅｓｓ．Ｔｈｅ ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｔｈｅ ｎｏｄｅｓ ｉｓ ｇｏｔｔｅｎ ｂｙ ａ ｖａｒｉａｔｉｏｎａｌ ｐｒｏｂｌｅｍ，ｗｈｏｓｅ Ｅｕｌｅｒ－Ｌａｇｒａｎｇｅ ｅｑｕａｔｉｏｎ ｉｓ ｔｈｅ ｎｅｗ ｍｏｖｉｎｇ ｍｅｓｈ ｅｑｕａｔｉｏｎ．Ｗｈｙ ａｎｄ ｈｏｗ ｔｏ ｆｏｒｍ ｔｈｅ ｎｅｗ ｍｅｔｈｏｄ ｉｓ ｉｎｔｒｏｄｕｃｅｄ ｉｎ ｄｅｔａｉｌｓ．Ｎｕｍｅｒｉｃａｌ ｒｅｓｕｌｔｓ ｉｎ ｔｗｏ—ｐｏｉｎｔ ｂｏｕｎｄａｒｙ ｐｒｏｂ－ ｌｅｍｓ ａｎｄ ｂｏｕｎｄａｒｙ ｌａｙｅｒ ｐｒｏｂｌｅｍｓ ａｒｅ ｖｅｒｙ ｅｎｃｏｕｒａｇｉｎｇ ａｎｄ ｃｌｅａｒｌｙ ｔｏ ｓｈｏｗ ｔｈａｔ ｉｎ ｍａｎｙ ｃａｓｅｓ ｔｈｅ ｍｅｔｈｏｄ ｒａｎ ｒｎｎｔｒｏ】ｔｈｅ ｍｅｓｈ ｍｏｔｉＯｎ． Ｋｅｙ ｗｏｒｄｓ：ｍｏｖｉｎｇ ｍｅｓｈ，ｖａｒｉａｔｉｏｎａｌ ｍｅｔｈｏｄ，ｍｏｎｉｔｏｒ ｆｕｎｃｔｉｏｎ，Ｅｕｌｅｒ—Ｌａｇｒａｎｇｅ ｅｑｕａｔｉｏｎ ０ 引言 面对大量复杂的非线性问题，因其中大部分问题的解在局部变化剧烈甚至出现间断，实际问题又迫切 需要尽可能精确掌握在此“奇异点”附近解的性态，所以在目前还没有一般方法解析求解非线性微分方程 的情况下，这就对数值解提出了极大的挑战．仔细分析这些非线性问题，可以发现它们与线性问题的最大 差异是大部分非线性问题的解只在局部有大梯度变化，而在其它大部分区域相对平滑．为了保证计算结果 的整体精度，一个自然的想法是在解变化剧烈的地方，网格分布得密一些，同时为了使计算量不至于增加 太多，在解变化比较平缓的地方，网格分布相对稀疏一些．对于发展型问题，“奇异点”一般随时间变化，网 格自然也应该随着时间运动以适应解的变化，即在解变化剧烈的地方保持细密的网格，这样既可以保证求 解的精度，又可避免付出过大的计算量，还可以增加时间步长．这就是运动网格的基本出发点． 运动网格方法已有３０年左右的研究历史，目前已得到很大的发展．无论对椭圆型问题、双曲型问题还 是抛物型问题，运动网格都有广泛的应用．运动网格可以与有限元方法相结合，也可与差分方法相结合．依 收稿日期：２００６－０９－２１．修回日期：２００６—１２—１３． 基金项目：国家９７３重点基础研究发展计划（２００５ＣＢ３２１７０Ｘ）资助项目． 作者简介：李征（１９７９一），女，硕士，主要从事计算数学方法的研究．Ｅ－ｍａｉｌ：ｌｉ—ｚｈｅｎｇ＠ｉａｐｃｍ．ａｃ．ｃｎ 通讯联系人：王双虎（１９６１一），研究员，主要从事计算数学的研究．Ｅ－ｍａｉｌ：ｗａｎｇｓｈｕａｎｇｈｕ＠ｉａｐｃｍ．ａｃ．ｃｎ 一１５— 维普资讯 http://www.cqvip.com 南京师大学报（自然科学版） 第３０卷第２期（２００７年） 据目前研究状况，按照所求解问题的不同，运动网格可分为两类：定态运动网格和动态运动网格． 运动网格最早的工作可以追溯到ｄｅ Ｂｏｏｒ¨ 的等分布原理．其主要想法是，如果网格的分布使解的某 种误差度量在计算网格上分布得比较均匀的话，那么采用这样的网格，将会在同等工作量的前提下获得较 好的结果．Ｗｈｉｔｅ 通过引进弧长坐标变换，将等分布原则用于计算两点边值问题，这是等分布原理最直接 的应用．Ｈｕａｎｇ 提出了运动网格偏微分方程（ＭＭＰＤＥ），将等分布原理应用到发展型问题，以ＭＭＰＤＥ控 制网格的分布．Ｍｉｌｌｅｒ 及其合作者将网格运动的思想与有限元相结合，提出了运动有限元方法（ＭＦＥ）， 并进一步发展为梯度加权运动有限元方法（ＧＷＭＦＥ） －ｓ］，而且给出了其变分解释和几何力学解释．也有 一类运动网格方法是在ＡＬＥ（Ａｒｂｉｔｒａｒｙ Ｌａｇｒａｎｇｅ Ｅｕｌｅｒ）框架下建立的，本质上可称之为拟Ｌａｇｒａｎｇｅ方法，即 利用Ｌａｇｒａｎｇｅ的思想进行网格运动控制．Ｃａｏ 提出了基于几何守恒律（Ｇｅｏｍｅｔｒｉｃ Ｃｏｎｓｅｒｖａｔｉｏｎ Ｌａｗ，ＧＣＬ） 的方法．由ＧＣＬ得到了网格速度场的散度条件，并根据Ｈｅｌｍｈｏｈｚ定理，添加适当的旋度条件来确定网格 速度场，与原方程联立可同时确定方程解及其网格． 本文得出了一种基于变分的运动网格方法．从网格应适应解的局部性态并具有一定光滑性的原则出 发，导出控制网格运动的泛函，并通过求解该泛函的极值问题确定网格的分布．在实际计算中，利用求解与 该泛函对应的欧拉一拉格朗日（Ｅｕｌｅｒ．ａｇｒＬａｎｇｅ）方程实现网格的重新分布．为了验证其正确性与有效性， 本方法应用于求解两点边值问题，以及难以用均匀网格计算的边界层问题，都得到了满意的数值结果． １网格变分问题的导出 微分方程数值解的本质就是数值积分，所以，不失一般性，考虑积分问题， ＝ ｂｕ（ ） 的数值计算． 首先，对区间［ａ，ｂ］进行网格剖分ａ＝ ＜ ＜…＜ ＜…＜ ＜ ＝ｂ，为了简单起见，选择矩形 积分公式Ｓ ＝∑ 一 （ 一 一 ）进行计算，其中， 一 ＝ｕ（ 一 ）． Ｊ＝１ 如图１所示，曲线为原函数，阴影部分表示数值积分 Ｓ ，空白部分的面积是数值误差．显然在确定数值积分公式 的情况下，网格的分布直接影响误差的大小．当选择均匀网 格计算时，得到的误差未必是最小的．我们希望在既不增加 计算结点数目，也不改变数值积分公式的前提下，通过合理 分布网格使计算得到的ｓ 与积分值， 尽可能相近，使误差 Ｅ１有 ｂ １｝， Ｅ ＝ｍｉ“…ｕ（ｘ）ｄｘ一∑ｕｊ一 （ 一 一 ）Ｉ Ｊ（１・１） 假定ｕ（ｘ）是解析函数，将其在每个区间［ ， ］上进 行Ｔａｙｌｏｒ展开： 均匀的网格上计算积分 优化的网格上计算积分 图１不同网格计算误差的比较 Ｆｉｇ．１ Ｎｕｍｅｒｉｃａｌ Ｅｒｒ ｆｏｒ ｄｉｆｅｒｅｎｔ ｍｅｓｈ ｆ ｕ（ｘ）ｄｘ一∑ｕｊ一 （ 一 一 ）＝了１∑［ｕ （ ． ）（ 一 ． ） ］＋Ｏ（Ａｘ ）． 可以看出，数值积分的误差主要来自于第一项．记作Ｒ ＝ １∑［ｕ （ 一 ）（ 一 一 ） ］． 引入计算参考区域［０，１］，作等间距剖分，０＝ ＜ ＜…＜ ＜…＜ 间距 ．（１．２） ＜ ＝１，其中 ＝ ， ＝ １由于 是常数，所以选择结点 使Ｒ 达到最小，等价于使Ｒ＝ Ｒ 达到最小，令 ＝ （ ），可以得到 Ｒ＝ 骞 （ Ｒ～ 、ｕ 蛙． ． 容易看出，当ｎ很大时， １６— 维普资讯 http://www.cqvip.com 李 征，等：一种基于变分的网格运动方法及其在边界层问题数值求解中的应用 进一步放大 ， ・・ ≤｝Ｊ＿ ［ｃ ＋（詈）‘］ ， 去掉常数，并不影响网格函数的确定．所以有理由认为，以下面变分问题的解所确定的网格，在其上采用矩 形积分公式（ ＝∑ｕｊ一 （ 一ｘｊ一 ））计算的积分值误差较小， ，＝ｍｉｎ『［（ｕｆ） ＋ （ ） ］ ． 其中，参数０≥０，可以根据实际问题的需要选取． （１．３） （１．３）式的意义很清楚：应该选择这样的网格，一方面在函数ｕ变化剧烈的地方要多分布网格，另一 方面又不应使网格的疏密变化过于剧烈．参数０起平衡的作用，用于控制网格变化的剧烈程度． 上面得到了在０为常数的情况下控制网格运动的变分形式（１．３），其实０也可根据问题需要不取为常 数．譬如， ７＝ｍｉｎ Ｊ［（ ） ＋（ （ ，ｕ）） ］ ． 对于常微分方程组： Ｙ １＝ （Ｙ１，Ｙ２，…，Ｙ ， ） （１．４） Ｙ ＝Ａ（ｙ ，ｙ ，…，ｙｍ，Ｘ） Ｙ ＝ｆｍ（Ｙ１，Ｙ２，…，Ｙ ， ）， 可以类似于单个方程的情况，将泛函取作 （１．５） ，＝Ｊ ｎ［（∑ （ …Ｙ…，Ｙ ，Ｙ ，…，）， ）（ ） ＋ ） ； ＝１ （１．６） ２ 网格控制方程 网格分布的确定问题已转化为求解泛函极小的问题．现在考虑如何求解此变分问题．目前有两种方 法，一种是直接求解（１．４），一种是求解其欧拉一拉格朗日方程．直接求解法的程序实现起来要复杂一些， 所以，采用差分方法求解其欧拉一拉格朗日方程． （１．４）对应的欧拉一拉格朗日方程为： ６１＝（ 一 ｄｏ）［（ ） ＋（ （ ）） ］＝一２［（ …（ ＋ ）） ２＋（Ｕｘ２＋０２） ］・ （２．１） 这就得到了控制网格质量的微分方程： （１．１，ｘｕ ＋０（０ ＋Ｏｕｕ ）） ；＋（ｕ ＋０ ） 描＝０． 如果原问题的求解区间为［ｎ，ｂ］，则方程（２．２）可以提边界条件： （０）＝ ， （１）＝ｂ． （２．２） 如果０为常数，则（２．２）可取为更简单的形式， １．１．“ ｘ；＋（ ＋ ） 描＝０． （２．３） 下面对方程（２．３）进行一些简单分析． （１）当Ｉ ｕ Ｉ很小时，网格控制方程的第一项可以忽略，（２．３）式近似为 硭＝０．考虑到边界条件，它 有惟一光滑解 ＝ｎ（１一 ）＋ ，确定的网格是均匀网格．其实不难想象，在Ｉ ｕ Ｉ很小时，均匀网格就是 合理的网格． （２）当０很大时，（２．３）式也近似为 ＝０，所以确定的网格也是均匀的．由此可见０控制了网格变化 的光滑程度，起正则化的作用． 为了协调（２．３）中第一项与第二项的关系，可以在方程（２．４）中引人参数Ａ，即 １．１，ｘｕ ；＋Ａ（ｕ：＋０） 描＝０，ｘ（Ｏ）＝０， （１）＝ｂ． 一（２．４） １ ７— 以后也会用（２．３）来控制网格． 维普资讯 http://www.cqvip.com 南京师大学报（自然科学版） 第３０卷第２期（２００７年） ３ 数值方法 下面考虑欧拉一拉格朗Ｅｔ方程（２．３）的数值求解．由于需求解的是非线性两点边值问题，数值求解 比较困难，这里给出两种求解方法． 第一种，根据（２．３）的特点，可以将 作为未知函数，采用坐标变换将非线性方程变成线性方程进行求 解．这样就得 “ “ 一（“ ＋０） ＝０． （３．１） 由于实际应用中空间网格时常是非均匀的，所以一阶导数和二阶导数的离散格式为： ： 霉 ： ， ２ ￣２ ［Ｌ （ （３．３） ＋ 一 ）（ ＋ 一 一 ）’（＋ 一 ＋ ）（ 方程（３．１）离散后，形成三对角系数矩阵是采用追赶法 求解，计算量不大． 下面需要利用得到的 ＝ （ ）确定 ＝ （ ），生成 新的物理空间网格如图２所示，求解方程（３．１）得到・ 对应的点． ： 采用抛物线来拟合函数．首先，选出与 ，＝上距离 ／／， 最近的 ｆ，其相邻的三个结点空间 、 ，和专＋ 可以惟 物理坐标 一确定所需的抛物线．对应于｝，值的元 即为所求的新物 图２确定新的物理空间网格示意图 理空间网格 ｝，即 对应的各点． Ｆｉｇ．２ Ｇｅｔｔｉｎｇ ｎｅｗ ｍｅｓｈ ｉｎ ｐｈｙｓｉｃａｌ ｄｏｍａｉｎ 另一种方法是时间相关法，就是把（２．３）的求解变为对非线性抛物型方程（３．４）的求解． “ ＝“ “ ；＋（“：＋０） ， （３．４） 初边值条件为： （ ，０）＝０＋（ｂ一０） ，ｘ（ｏ，ｔ）＝０， （１，ｔ）＝ｂ． （３．５） 可以认为当ｔ一∞，（３．４）求得的解趋于（２．３）的解．由于网格的好坏只是一个相对概念，所以对于 实际问题而言，只要求解到ｔ足够大即可，不一定非要达到收敛． ４ 数值试验 为了验证网格运动方程的有效性，以常微分方程的两点边值问题 Ｉ（ ）＝ ，“ （４．１） ．ｕ（０）＝Ｏｔ， “（ｂ）＝ 作为算例进行检验．这里给出了两个典型算例的计算结果．采用打靶法对两点边值问题进行求解，每迭代 一步，先数值求解两点边值问题，再依照数值解生成新的网格，循环求解，直到收敛． Ｍｏｄｅｌ １ “”＋２ｕｕ ＝０，ｌＬ（０）＝１，ｌＬ（１）＝２． （４．２） 将其写成方程组的形式： ，２ ．， （４．３） 方程的精确解“（ ）＝ Ｋ－／，，￣ｅｘｐ ＺＯｔｅｘｐ（ｌ２ａｘ／Ｘ／ ｅ）１（ ，、， Ｊ＋上）＋１ 一，其中，Ｋ＝ 一上 ！ 二 ２ １ ± ２ Ｏｔ ， 满足ｅｘｐ（＼ ）＝ ， （Ｏｔ＋１）（Ｏｔ一２） 维普资讯 http://www.cqvip.com 李 征，等：一种基于变分的网格运动方法及其在边界层问题数值求解中的应用 （ ＋ ） ；＋Ｔ（ｆ ＋厂；＋ｏ）ｘ转＝０， （４．４） 其中， ＝ｖ／８， ＝一２ｕｖ／８， ＝一２ｕｖ／８ ， ＝一２（ｆ￣ｖ＋Ｌｕ）／ｅ ． 求解区间离散为１０个单元，共有 １１个结点．图３给出了 ＝１， ：１时 的数值结果．（ａ）给出了均匀网格和 自适应网格上数值解的误差比较图． 其中，方形表示均匀网格上数值解的 误差，圆形、菱形、三角为采用我们的 运动网格方法在０取不同值时数值解 的误差．（ｂ）为０：１时数值解与网 格关系图．其中，实线表示精确解，方 （ａ）不同网格计算的误差比较图 （ｂ）０＝１时数值解与网格分布图 块表示数值解，纵线标出了网格结点 图３算例１中ｅ＝ｌ。３，－－１的数值结果 Ｆｉｇ．３ Ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ｗｉｔｈ ｅ＝ｌ，ｙ ｌ ｏｆ ｍｏｄｅｌ １ 的位置，可见，网格结点在解变化剧烈 的地方汇聚． 当 减小的时候，本问题会在 ＝０附近出现边界层效应， 越小，边界层效应越明显．如果采用均匀网 格计算这类问题，除非网格非常细密，否则不能描绘出边界层内部的信息．图４给出了对应于 取不同值 的网格与数值解．实线代表精确解，圆点的横坐标为网格结点的位置，纵坐标代表数值解．参数 ＝２，０＝ ０．５．其中，图４（ａ）给出了 分别取１０～、１０～、１０ 和１０ 的数值解和网格分布关系图．随着 一０，边界 层越来越窄，网格大部分汇聚在边界层内，图４（ｂ）给出了对应于图４（ａ）的边界层内部的网格分布，相当 于各图的局部放大．表１给出取不同占值的最大模误差和平方误差．网格的分布使误差没随着边界层效应 ≈ 的剧烈而增加．这个例子说明，；｝｝ ２ ９ ８ ７ ６ ５ ４ ３ ２本方法可以用于求解此类问题． １ ｌ ２ １．８ １．６ １．４ １．２ １ ０ ０，２ ０．４ ０，６ ０．８ １ ０ ０．２ ０．４ ０．６ ０．８ １ ０ ０．２ ０．４ ０．６ ０．８ １ ０ ０．２ ０．４ ０．６ ０．８ １ ｘ ｘ ｘ ｘ （ａ）不同ｓ的数值解和网格图 （ｂ）边界层内部的网格分布图 图４算例１中 取不同值的数值结果 Ｆｉｇ．４ Ｎｕｍｅｒｉｃａｌ ｒｅｓｕｌｔｓ ｏｆｄｉｆｆｅｒｅｎｔ ８ 表１ ８取不同值的误差比较 ＪＴａｂｌｅ １ Ｎｕｍｅｒｉｃａｌ ｅｒｒ ｏｆ ｄｉｆｆｅｒｅｎｔ ｓ 【． 罟ｄ ＝ｖ＝ ，ｕ １ ｃ。 ，． ，ｕ（１）＝ｌ ４．５ １０一。４ ５２５４４９１０５５２４９２６Ｅ一００２ ２．８７７５８３６６６６１ １５３３Ｅ—ｏｏ２ １０— ４．５４３４２０３７９７４９５７４Ｅ一００２ ２．８８６６２１７３４４７６９７８Ｅ一００２ １０一 １．９０７４７０９４５５４６５５２Ｅ一００２ ｌ＿１３４０６０３５４６７３４７５Ｅ一００２ １０—８ ３．８８８１９７４０５０２１５１０Ｅ一００２ １．９１５４０８１９３９９３９５６Ｅ一００２ 维普资讯 http://www.cqvip.com 南京师大学报（自然科学版） 第３０卷第２期（２００７年） 运动网格方程： （／Ｊｌ ＋ ） ；＋ＴＣ＋ ＋０）ｘ菇＝Ｏ， 其中， ＝ ， ＝ ， ＝一 １，２ ＝一 １，． （４．６） １ ０．９ ０．８ ０．７ ０．６ ０．５ ０．４ ０－３ ０．２ ０．１ 图５ ｅ－－０．５时 取不同值的数值结果 Ｆｉｇ．５ Ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ｏｆ ｄｉｆｆｅｒｅｎｔ ０ ｗｈｅｎ ｅ＝０．５ 图５给出了 ＝Ｏ．５， ＝１时的数值结果，图５（ａ）为均匀网格与运动网格方法得到的数值解的误差 比较图．（ｂ）为０＝１时数值解与网格关系图．同样可见，自适应网格的数值误差要小于均匀网格． ０ ０ ￡０ ０ 图６不同￡的运动网格和均匀网格计算的误差比较图 Ｆｉｇ．６ Ｎｕｍｅｒｉｃａｌ ｅｒｒ ｏｆ ｍｏｖｉｎｇ ｍｅｓｈ ａｎｄ ｕｎｉｆｏｒｍ ｍｅｓｈ ｆｏｒ ｄｉｆｅｒｅｎｔ￡ 当 减小的时候，解在 ＝Ｏ处出现边界层．当 很小的时候，１ １个结点已经不足以给出适当的解，因 此，随着 的改变选择不同的结点个数求解．图６给出了 取不同值时，运动网格和均匀网格的误差比较 图．圆点代表的运动网格误差要远远小于均匀网格的误差．表２给出了在取不同 值Ⅳ个结点情况下的最 （ａ） 取不同值时方程解与网格的关系图 ｆｂ ８＝０．ｂ１０１和８－－０．００１时边界层内部的网格分布图 图７￡取不同值的数值结果 Ｆｉｇ．７ Ｎｕｍｅｒｉａｌｃ ｒｅｓｕｌｔｓ ｏｆｄｉｆｅｒｅｎｔ 一２０一 维普资讯 http://www.cqvip.com 李 征，等：一种基于变分的网格运动方法及其在边界层问题数值求解中的应用 大模误差和平方误差图７给出了 分别取０．２、０．１、０．０１和０．００１的计算结果．（ａ）给出了 取不同值时 网格与解的关系图，（ｂ）给出了 ＝０．０１和 ＝０．００１时候边界层内部网格与解的关系图．参数ｙ＝０．５， ０＝０．１． 表２ ｓ取不同值的误差比较 Ｔａｂｌｅ ２ Ｎｒｉｍｅｒｉｃａｌ ｅｒｒ ｏｆ ｄｉｉｅｒｅｎｔｆ ｓ ５结论 本文给出的基于变分形式的运动网格方法可以自然推广到高阶方程、方程组，以及偏微分方程．作为 对该运动网格方法的检验，应用此方法于一系列有不同边界层效应的两点边值问题的数值模拟，结果说明 该方法给出了合适的网格分布． ［参考文献］ ［１］ Ｃ ｄｅ Ｂｏｏｒ．Ｇｏｏｄ ａｐｐｒｏｘｉｍａｔｉｏｎ ｂｙ ｓｐｌｉｎｅｓ ｗｉｔｈ ｖａｒｉａｂｌｅ ｋｎｏｔｓ ＩＩ［Ｍ］／／Ｓｐｒｉｎｇｅｒ Ｌｅｃｔｕｒｅ Ｎｏｔｅｓ Ｓｅｒｉｅｓ ３６３．Ｂｅｄｉｎ：Ｓｐｒｉｎｇｅｒ— Ｖｅｄａｇ，１９７３． ［２］ Ｗｈｉｔｅ Ａ Ｂ．Ｏｎ ｓｅｌｅｃｔｉｏｎ ｏｆ ｅｑｕｉｄｉｓｔｒｉｂｕｔｉｎｇ ｍｅｓｈｅｓ ｆｏｒ ｔｗｏ—ｐｏｉｎｔ ｂｏｕｎｄａｒｙ—ｖａｌｕｅ ｐｒｏｂｌｅｍｓ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｎｕｍｅｒｉｃａｌ Ａｎａｌｙｓｉｓ，１９７９，１６（３）：４７２—５０２． ［３］Ｈｕａｎｇ Ｗ，Ｒｅｎ Ｙ，Ｒｕｓｓｅｌｌ Ｒ Ｄ．Ｍｏｖｉｎｇ ｍｅｓｈ ｐａｒｔｉａｌ ｄｉｆｆｅｒｅｎｔｉｌａ ｅｑｕａｔｉｏｎｓ（ＭＭＰＤＥｓ）ｂａｓｅｄ ｏｎ ｔｈｅ ｅｑｕｉｄｉｓｔｒｉｂｕｔｉｏｎ ｐｒｉｎｃｉ— ｐｉｅ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｎｕｍｅｒｉｃａｌ Ａｎａｌｙｓｉｓ，１９９４，３１：７０９．７３０． ［４］Ｍｉｌｌｅｒ Ｋ，Ｍｉｌｌｅｒ Ｒ Ｎ．Ｍｏｖｉｎｇ ｉｆｎｉｔｅ ｅｌｅｍｅｎｔ Ｉ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｎｕｍｅｒｉｃａｌ Ａｎａｌｙｓｉｓ。１９８１，１８：１ ０１９．１ ０３２． ［５］Ｍｉｌｌｅｒ Ｋ．Ｍｏｖｉｎｇ ｉｆｎｉｔｅ ｅｌｅｍｅｎｔ ＩＩ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎｌａ ｏｎ Ｎｕｍｅｒｉｃａｌ Ａｎａｌｙｓｉｓ，１９８１，１８：１ ０３３—１ ０５７． ［６］ Ｃａｒｌｓｏｎ Ｎ，Ｍｉｌｌｅｒ Ｋ．Ｄｅｓｉｇｎ ａｎｄ ａｐｐｌｉｃａｔｉｏｎ ｏｆ ａ ｇｒａｄｉｅｎｔ—ｗｅｉｇｈｔｅｄ ｍｏｖｉｎｇ ｉｆｎｉｔｅ ｃｏｄｅ，Ｐａｒｔ Ｉ，１一Ｄ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｓｃｉｅｎｔｉｆｉｃ Ｃｏｍｐｕｔｉｎｇ，１９９８，１９：７２８—７６５． ［７］Ｃａｒｌｓｏｎ Ｎ，Ｍｉｌｌｅｒ Ｋ．Ｄｅｓｉｇｎ ａｎｄ ａｐｐｌｉｃａｔｉｏｎ ｏｆ ｇｒａｄｉｅｎｔ—ｗｅｉｇｈｔｅｄ ｍｏｖｉｎｇ ｆｉｎｉｔｅ ｃｏｄｅ，Ｐａｒｔ ＩＩ，２一Ｄ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｓｃｉｅｎｔｉｆｉｃ Ｃｏｍｐｕｔｉｎｇ，１９９８，１９：７６６－７９８． ［８］ Ｍｉｌｌｅｒ Ｋ．Ａ ｇｅｏｍｅｔｉｒｃａｌ—ｍｅｃｈａｎｉｃａｌ ｉｎｔｅｒｐｏｌａｔｉｏｎ ｆ ｏｒａｄｉｇｅｎｔ—ｗｅｉｇｈｔｅｄ ｍｏｖｉｎｇ ｉｆｎｉｔｅ ｅｌｅｍｅｎｔｓ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｎｕｍｅｒｉ— ｃａｌ Ａｎａｌｙｓｉｓ，１９９７，３４：６７－９０． ［９］Ｃａｏ Ｗ，Ｈｕａｎｇ Ｗ，Ｒｕｓｓｅｌｌ Ｒ Ｄ．Ａ ｍｏｖｉｎｇ ｍｅｓｈ ｍｅｔｈｏｄ ｂａｓｅｄ ｏｎ ｔｈｅ ｇｅｏｍｅｔｉｒｃ ｃｏｎｓｅｒｖａｔｉｏｎ ｌａｗ［Ｊ］．ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｓｃｉ— ｅｎｔｉｌｔｅ Ｃｏｍｐｕｔｉｎｇ，２００２，２４（１）：１１８—１４２． ［责任编辑：陆炳新］ 一２１—