非均匀三角形网格下椭圆问题的代数两网格法

高校应用数学学报　２０１３，２８（４）：４３９—４４６　

非均匀三角形网格下椭圆问题的　

代数两网格法　

李　明　，李郴良。　

ｆ１．红河学院数学学院，云南蒙自６６１１００；　

２．桂林电子科技大学数学与计算科学学院，广西桂林５４１００４）　

摘要：针对非均匀三角形网格离散带各向异性系数或间断系数的二维椭圆问题形成　

的多尺度系统，基于极大不相关子集的粗化算法，构造一种插值算子，提出了一种代数　两网格法．数值实验表明新算法的有效性．　

关键词：非均匀三角形网格；椭圆问题；插值算子；代数两网格法　

中图分类号：Ｏ２４１．６　

文献标识码：Ａ　文章编号：１０００．４４２４（２０１３）０４—０４３９—００８　

§１　引　言　

多重网格法是求解偏微分方程常用的一类数值方法，可分为几何多重网格（ＧＭＧ）法和代　数多重网格（ＡＭＧ）法．相对于ＧＭＧ法，ＡＭＧ法几乎不需要使用问题的几何信息，能够有效处　

理很多复杂问题，包括非结构网格、各向异性系数和变系数等问题，但对于强震荡、间断系数　

和自适应、大变形网格等问题，经典的ＡＭＧ法（见『１—３］）难以有效求解．文献ｆ４—５］指出上述不足　

的原因是这些问题导出的矩阵一般具有如下性质：同一行非对角元素在数值上相差很大　甚至　相差几个、几十个数量级，并称具有这种性质的线性系统为多尺度系统ｆ否则称为单尺度系统）．　

针对这一情况，徐小文等人通过将系数矩阵中某些特殊的单尺度性质的连通子块找出，形成单　

尺度子块，然后基于这些连通子块构造特殊的粗化和松弛策略，提出一种基于局部松弛和粗化　

的ＡＭＧ法（见ｆ４—５１）．谭敏等人利用网格点之间的几何距离，将所有点分离为各向异性部分和各　

向同性部分，然后再进行粗化，提出了一种各向异性四边形网格下的代数多重网格法（见ｆ６］）．但　

此方法需要使用网格节点的几何坐标信息计算节点间的距离．　

为了进一步发展不依赖于多尺度性质的ＡＭＧ法，本文针对非均匀三角形网格离散二维椭圆　

问题形成的有限元方程，使用极大不相关子集（ＭＩＳ）来生成粗点集，构造相对应的插值算子，提　

出一种代数两网格法．新算法无需形成单尺度子块，也不用讨论各向同性或各向异性以及节点　

收稿日期：２ｏ１２—０８—２０　修回日期：２０１３—０５—０９　基金项目：国家自然科学基金（１１１６１０１４）；云南省科技厅青年项目（２０１２ＦＤ０５４）；红河学院硕博项目（ＸＪ１Ｓ０９２５）

　高校应用数学学报　第２８卷第４期　

间的几何距离等情况，新算法流程更为简单．数值实验表明新算法在细层上无需磨光步，具有计　

算量少，计算时间短，稳健性强的优点．　

考虑如下二维椭圆问题　

一　笔　～　＝，，　在　，　ｌ　：０，　在ａ　，　（１）　

其中，是己知足够光滑的函数，　，　为大于零的己知系数，ｒ“为足够光滑待求函数，　是一个二　

维凸多边形区域，ａ　是　的边界．　

使用非均匀三角形网格剖分区域　，采用线性拉格朗曰有限元进行离散问题（１），记对应的　

有限元方程为　

Ａ１Ｕ１＝Ｆ１　

其中Ａｌ＝（。巧）　是咒阶对称正定方阵，　＝（＾，ｆ２，・　向量．　

§２代数两网格法　

粗化算法和插值算子是ＡＭＧ法的重要组成因素．粗化算法的粗化效果对ＡＭＧ法的计算效　

率有着重要的影响作用．为了快速形成粗点集，本文采用极大不相关子集（ＭＩＳ）粗化算法．具体　

如下　

定义２．１［７］在网格图Ｇ１中，如果某一（顶点）指标集Ｃ１满足如下两个条件，则称之为Ｇ１的极　

大不相关子集（ＭＩＳ）：　

条件１．对　中任意一对指标，相应的顶点之间没有边相连．　

条件２．当增加任意一个指标到Ｃ１中时，　１不再满足条件１．　

记有限元方程（２）对应的变量指标集合口＝｛１，２，・一，几），使用图论的知识，根据矩阵　１可　以构造与　１对应的网格图Ｇ１．由定义２．１可得Ｇ１的极大不相关子集（ＭＩＳ），将其记为粗网格变量　

指标集Ｃ１，并将　１中元素按升序排列．记口关于Ｃ１的补集为Ｓ１＝ｎ＼ｃ１．　

为了便于下面讨论，引入如下标记　

＝｛Ｊ：１ａｉｊｌ≠０，ｉ≠　）为ｉ点的邻点构成的集合．　

孵＝Ｃ１　ｎ　为ｉ点的邻点中粗点构成的集合．　

＝Ｓ１　ｎ　为ｉ点邻点中细点构成的集合．　

ｓ＝ｉｎｄ（ｉ）为粗点指标ｉ∈Ｃ１在粗点集Ｃ１中的序号．　

Ｉｃ１　Ｉ表示粗点集Ｃ１中元素个数．　

１　Ｉ表示集合　中元素个数．　

易知，集合Ｃ１，Ｓ１，　非空．由定义２．１可知，若ｉ∈　１，则Ⅳ　非空．事实上假设　为空集，　

Ｎｉ点的所有邻点中都没有粗点，这意味着ｉ作为粗点时，也不会破坏定义２．１条件１，由定　义２．１知　∈Ｃ１，这与ｉ∈Ｓ１矛盾．　

下面讨论构造基于ＭＩＳ的插值算子．方程（２）可写为　

∑。　＾，　

Ｊ＝１　、　Ｕ　２歹　维　佗　是　Ｔ　、ｌ，　ｎ　２　Ｌ　／●　＝　ｌ　Ｔ　、Ｉ，　

Ａ　李　明等：非均匀三角形网格下椭圆问题的代数两网格法　４４１　

即　ａｉｉｕｉ＋∑ａｉｄｕｊ＝＾，ｉ＝ｌ　２一，ｎ　

ｊＣＮｉ　假设粗点集　对应的解　＝［ＵＨ，１，　日，２，　

之对应点的近似值．即　ＵＨ，ＩＣ　ｆ］Ｔ已经求出，则将其作为细网格上与　

７２ｉ：≈　日．ｉｎｄ（０，Ｖｉ∈Ｃ１　

因　＝　Ｕ　，对于Ｖｉ∈ｓ１，式（３）可写为　

ｎｉｉ？￣ｉ＋　二ｎｉｊｕｊ＋　二０“　ｔ＝　，Ｖｉ∈ｓ１　

ｊｅＮ．：ｔＥＮ　２　

因为＾　非空，对上式左端∑ａｉｒｕｔｄ￣的　ｔ做如下近似　ｔ∈Ⅳ　

再将式（２）和式（４）代入式（３），移项得　

由４１正定可知，　∑　

≈　，ｖｔ∈　≈商’吼　

％　ｔ≈＾一∑吼ｊＵＨ　）　

ｊ∈Ｎ；　ｔ　ＵＨ，ｉｎｄ（ｊ）　

一　，　∈　Ｉ　Ｉ　…　

（一　Ｚ　ａｉｊｕＨ，ｉｎｄ（ｊ）ｏＮ？）ｊ　Ｖｉ　Ｅ￥１　一　，

将式（２）和式（５）综合起来　

式（６）可记为　

其中　札日，ｉｎｄ（ｉ）　

上ａｉｉ（一　Ｅ∑Ｎｃ　＋　＋　

ｆ　

（Ｐｉｊ）似ＩＣ　ｆ＝｛　【　

＂ｆｌｕ１＝　ｕｏ：≈　＋Ｗ７１　

ａｉｋＩ．Ⅳ　ｌ＋∑ａｉｔ　

Ⅱ。ｔｌⅣ　ｌ　’　

ｃ　＝　姜　

基于上述粗化算法和插值算子，结合文献［８－ｌＯ］的思想，　）　若ｉ∈Ｃ１　若ｉ∈Ｓｌ　

若ｉ∈Ｃ１，Ｊ＝ｉｎｄ（ｉ）　

若ｉ∈Ｓ１，ｋ∈Ⅳ　，Ｊ＝ｉｎｄ（ｋ）　

其它　

若ｉ∈Ｃ１　

若ｉ∈Ｓ１　

提出如下算法　（２）　

（５）　

（６）　％　

≈　４４２　高校应用数学学报　第２８卷第４期　

图ｌ　ｔｐｙｅｌ（左图，纵轴、横轴方向不同尺度）；ｔｐｙｅ２（￣Ｆ图，横轴方向不同尺度）；ｔｐｙｅ３（右图，纵轴方向不同尺度）　

算法１　代数两网格法　

步骤１．生成粗网格矩阵和右端量：ＡＨ：＝（　）ＴＡ１　，　：＝（　）Ｔ　Ｆｌ　

步骤２．直接求解ＡＨＵＨ＝ＦＨ得　；　

步骤３．插值：ｕ１０：＝Ｐ　＋　１；　

步骤４．若　＜￡１，令　：＝　，否则磨光　次：　：＝Ｇ　（ｕｏ）．　

其中Ｇ为磨光算子，本文取￡１＝１０＿。。．　

§３数值实验　

为了验证算法的有效性，本文选用的通常的代数二重网格法作为对比算法．其粗化算法与　本文粗化算法相同，其插值算子￣ｌｆｌＴ［１ｌ１　

Ｉ　１，　若ｉ　Ｅ　Ｃｉ，Ｊ＝ｉｎｄ（ｉ）　Ｐ＝（Ｐｉｊ）　ＩＧ　ｌ＝｛丽１，　若ｉ∈Ｓ１，ｋ∈　，Ｊ＝ｉｎｄ（ｋ），ｍ＝ｌ　ｌ　【

０，　其它　对比算法的前后磨光步数均取３，初始值为零向量，对比算法的终止条件为相邻两次近似　

解的误差小于１０－６或调用次数达￣，ＪｌＯ０次．本文算法与对比算法的磨光算子Ｇ为共轭梯度法．使　

用ＭＡＴＬＡＢ２００８软件编写数值实验的程序代码，并在处理器为２．４Ｇ，内存为２Ｇ的台式电脑上进　

行计算．　

采用用图ｌ中网格离散问题区域，使用三角形有限元离散算例１和算例２．约定如下标记：　“ｔｙｐｅ”为网格类型，图２和图３中，“Ａｌｇｏｒｉｔｈｍ　１”与“Ｃｏｍｐａｒｅ　ａｌｇｏｒｉｔｈｍ”分别表示本文算法１和　

对比算法．“　”表示非均匀网格中最大单元与最小单元的直径之比，“　，　”为问题（１）中的　

系数，ｎＵｍｇｒｉｄ为网格单元总数，ＩＩ札一　１Ｉｏ。为问题解札与算法求出的数值解　的无穷范数误差，　

“ｕ”为对比算法的迭代次数，“　”为算法１步骤４中磨光步数，“ｃｐｕ”为算法的计算时间（秒）．　

算例１各向异性问题　

右端量ｆ＝一　（一ｓｉｎ（ｙ）ｅ　（１一Ｘ２）（１一Ｙ。）＋４ｘ　ｓｉｎ（ｙ）ｅ　（１一Ｙ　）一２ｓｉｎ（ｙ）（１一ｅｘ）（１一ｙ２））一　

（一ｓｉｎ（ｙ）（１一ｅ　）（１－ｘ　）（１－ｙ　）－－４ｃｏｓ（ｙ）（１－－ｅ　）（１一　。）ｙ－２ｓｉｎ（　）（１一ｅ　）（１一　）），　≠　是　各向异性系数，真解　＝ｓｉｎ（ｙ）（１一ｅｘ）（１一Ｘ２）（１一ｙ２），区域　：（０，１）×（０，１）．　

算例２间断系数问题　

右端量，＝　丌２（１一ｃ０ｓ　（７ｒＸ）一ｓｉｎ（７ｒＸ）＿ｅＳｉｎ（Ｔｒｘ），。ｓ—ｉｎ（ｃＴ　ｒｙ）。，间断系数为　＝　＝｛　：　

真解钆＝ｓｉｎ（丌　）（南一１），区域　：（０，１）×（０，１）．

　李　明等：非均匀三角形网格下椭圆问题的代数两网格法４４３　

图２例题ｌ，ｔｐｙｅｌ￣，鲁　＝１００，Ｑ　＝１；左图：迭代步数口ｓ参数ａ　（ｎｕｍｇ　ｉｄ＝８１９２）；右图：迭代步数＂ｓ残　量（ａ　＝１０　，ｈｕｍｇｒｉｄ＝３２７６８）　

图３例题１，ｔｐｙｅｌ￣，Ｑ　＝１，　ｈｍ　ｉｎ＝１００，左图：ＩＩｕ—ｕ；ｌｌｏｏ　ｕｓ参数　（ｎｕｍｇｒｉｄ＝３２７６８）；中图：Ｉ｝ｕ—Ｕ｛ｌｌｏｏ　Ｖ８单元数ｎｕ　ｇｒｉｄ（　＝１０　）：右图：迭代步数＂ｓ单元数ｎｕｍｇｒｉｄ（ｏｒ　＝１０。）　

图２（左图）中，横轴表示参数　的取值，纵轴为算法１的磨光步数或对比算法的迭代次数．可　以看出，当参数　依次取１，１０，１０　，…，ｌ０　时，对比算法至少需要１００次迭代，才能得到预定精度　

的数值解．而本文算法１，不需要任何磨光步．这说明算法１的磨光步数不受参数　取值的影响．　

图２（右图）中，横轴表示迭代次数（或磨光步数），纵轴为表示算法１或对比算法对应的残　量．可以看出，对比算法消去残量的作用并不明显，而算法１无需磨光步，对应的残量近似为　

零ｆ≈１０－１０）．这说明算法１只需要插值就可以得到有限元方程的高精度近似解．　

图３（左图）中，横轴表示参数　的取值，纵轴表示算法１或对比算法求出的数值解　与问题　真解　的无穷范数误差．可以看出，当参数　依次取１０～，１０＿。，…，１０　时算法１的计算精度明显　

高于对比算法．　

图３（中图）可以看出，当网格单元数量增加时，即加密网格时，算法１的计算精度逐渐增加，明　

显优于对比算法．　

由图３（右图）可以看出，算法１的磨光步数不受求解规模的影响，明显少于对比算法的迭代次　

数．　

从表１至表４可以看出：　

１）磨光步数或迭代次数方面：　

算法１在细层上无需磨光步，也就是仅通过插值算子，将粗层上的解插值到细层，就可以得　到满足业　＜Ｅ１的数值解　．而对比算法的迭代次数至少需要１００次才能达到预定的计　

算精度．