椭圆方程数值解

j.椭圆方程数值解法本章考虑椭圆微分方程数值解法。首先以二维二阶椭圆方程为例，给出矩形网和三角网上的差分法。然后以一维二阶椭圆方程为例，简要描述有限元法的基本思想。J.1矩形网上差分方程考虑二维区域（区域=连通的开集）G上的二阶椭圆型偏微分方程第一边值问题（j.1）,,,xxyyxyuuCuDuEuFxyuxyxyG其中C,ED,是常数；0E；GC0,yxFF；(,)xy是给定的光滑函数；是G的边界；GG。假设（J.1）存在光滑的唯一解。考虑一种简单情形，即求解区域G是矩形区域，并且其四个边与相应坐标轴平行。令1h和2h分别为x和y方向的步长，用平行于坐标轴的直线段分割区域G，构造矩形网格：hG为网格内点节点集合，h为网格边界节点集合，hGhGh。对于内点jiyx,hG，用如下的差分方程逼近微分方程（J.1）：(J.2)1,1,,1,11,1,,1,122212122222ijijijijijijijijijijijijuuuuuuuuuuCDEuFhhhh其中),(jiijyxFF。（J.2）通常称为五点差分格式。方程（J.2）可以整理改写为（J.3）jia,1jiu,1+jia,1jiu,1+1,jia1,jiu+1,jia1,jiu+jia,jiu,ijF对每一内点jiyx,都可以列出这样一个方程。方程中遇到边界点时，注意到边界点上函数值u已知，将相应的项挪到右端去。最后得到以u的内点近似值为未知数的线性方程组。这个方程组是稀疏的，并且当1h和2h足够小时是对角占优的。用（J.1）的真解(,)uxy在网点上的值(,)ijuxy、1(,)ijuxy等等分别替换（J.2）中的iju、1,iju等等，然后在(,)ijxy点处作Tailor展开，便知差分方程（J.2）逼近微分方程（J.1）的截断误差阶为2221hhO。另外可以证明，五点差分格式的收敛阶为2212()Ohh，并且关于右端和初值都是稳定的。矩形网格差分格式的优点是计算公式简单直观。但是，当G是非矩形区域，并且边界条件包含法向导数（第二和第三边值条件）时，在矩形网格边界点建立差分方程是一件颇为令人烦恼的事情。矩形网格的另一个大缺点是不能局部加密网格。图J.1一般区域G的矩形网格J.2.三角网差分格式本节我们将积分插值法用于三角网，建立三角网差分格式。三角网差分格式具有网格灵活和法向导数边界条件易于处理等优点，特别地，它还保持积分守恒（质量守恒），深受使用者欢迎。文献上常称之为有限体积法或广义差分法。考虑有界区域G上的Poisson方程（J.4）fu，(,)xyG在边界的各个部分1、2和3分别给定第一、第二和第三边值条件：Oxy（J.5a）1(,)uxy（J.5b）),(yxu2Γn（J.5c）),(yxuu3Γn其中是常数，n是边界的外法向。作G的三角剖分：在上取一系列点，连成闭折线~，并记G~为由~围成且逼近G的多边形区域。将G~分割成有限个三角形之和，使每个三角形的每个内角不大于90，并且每个三角形的任一顶点与其他三角形或者不相交，或者相交于顶点。引入如下术语。节点：三角形的顶点；单元：每个三角形；相邻节点：同一条边上的两个节点；相邻单元：有一条公共边的两个三角形。对于任一节点，考虑所有以它为顶点的三角形单元和以它为顶点的三角形边，过每一条边作中垂线，交于外心，得到围绕该节点的小多边形，称为对偶单元。全体对偶单元构成区域G的一个新的网格剖分，称为对偶剖分。图J.2三角网及其对偶剖分0P~1P2P3P4P5P6P0P1q2q3q4q5q6q1m2m3m4m5m6m0G(a)0P1P2P3P4P1m2m3m4m1q2q3q0G(b)图J.3内点(a)与边界点(b)的对偶单元11mqusdn01011uuppqm1，12qqusdn02021uuppqq2，23qqusdn03032uuppqq3，34qmusdn340044qmuupp，对于supmdn04和sumpdn10，分别利用右矩形公式和梯形公式计算所涉及到的积分，导出如下差分近似：4404040404040002mmpmpmpmpmpusussusmpuuddddn440004000041222mmpuumpuuu40004324ppuu1010mpmpususddn10001324ppuu这里00()p。将上述六个公式带入（J.6）中，就得到边界点0p的差分方程。所有内点和边界点的差分方程构成一个封闭的线性方程组，其系数矩阵是稀疏的，并且当0时是对称的。J.3椭圆方程的有限元法有限元法是与差分法并驾齐驱的一套求解偏微分方程的方法。它的基本想法是，首先把微分方程转化成一种变分方程（微分积分方程），从而降低了对解的光滑性和边值条件的要求；然后，把求解区域划分成有限个单元（有限元），构造分片光滑函数，这个光滑函数由其在单元顶点上的函数值决定；最后，把这个分片光滑函数带入到上述微分积分方程中去，就得到关于单元顶点函数值的一个线性方程组，解之即得有限元解。与差分法相比，有限元法易于处理边界条件，易于利用分片高次多项式等等来提高逼近精度。函数集合mH作为例子，我们将考虑区间0,1I上的椭圆微分方程。用2()LI表示在I上勒贝格平方可积函数的集合，()mHI表示本身以及直到m阶的导数都属于2()LI的函数的集合。我们下面用到的主要是1()HI。这里所说的导数准确地说是应该是广义导数，对此我们不予详细说明，只需知道比如说，连续的分片线性函数（折线函数）就属于1()HI，其广义导数是分片常数函数。另外，我们还用到函数集合11(){(),(0)0}EHIvvHIv。变分方程考虑两点边值问题(),(0,1)puqufx（J.7a）(0)0u（J.7b）0)1(u（J.7c）其中,,pqf都是区间(0,1)上的光滑函数，并且0pp，0p是一个正常数。用1()EHI中任一函数v乘（J.7a）式两端，并在[0,1]上积分，得10[()]0puvquvfvdx（J.8）利用分部积分，并注意0)1(u和0)0(v，得101010|)(dxvupvupdxup10dxvup以此代入到（J.8）得到100)(dvfvquvvup（J.9）为了方便，定义10,wvwvdx（J.10）),(),(),(vqwvwpvwa（J.11）则相应于微分方程（1）--（3）的变分方程为：求1()EuHI满足),(),(vfvua)(1IHvE（J.12）注意在（J.12）中不出现二阶导数。我们已经看到，满足微分方程（J.7）的光滑解一定满足变分方程（J.12）。而变分方程（J.12）的解称为微分方程（J.7）的广义解，它可能只有一阶导数，因此可能不是（1）-（3）的解；但是如果它在通常意义下二阶可微，则一定也是（J.7）的解。另外，注意在变分方程（J.12）中，强制要求广义解u满足边值条件(0)0u，因而称之为强制（或本质）边界条件；而对边值条件0)1(u，则不加要求。但是可以证明，如果广义解u在通常意义下二阶可微，则一定有0)1(u，即这个边界条件自然满足。这类边界条件称之为自然边界条件。总之，变分方程（J.12）不但降低了对解的光滑性的要求，也降低了对边值条件的要求。有限元空间构造有限元法的第一步与差分法一样，也是对求解区间作网格剖分0101nxxx。相邻节点1,iixx之间的小区间1,iiiIxx称为第i个单元，其长度为1iiihxx。记maxihh。顺便说一下，有限元法不要求步长ih是常数。而差分法通常要求步长ih是常数，以免截断误差阶数降低。在空间1()EHI中，按如下原则选取有限元空间hV：它的元素()hux在每一单元上是m次多项式，并且在每个节点上都是连续的。当1m时，就得到最简单的线性元，这时每个hhuV可表为11(),iihiiiiixxxxuxuuxIhh，1,2,,in（J.13）其中(),ihiuux0(0)0huu。图J.3.一维线性元线性元的另外一种表示方法用到以下具有局部支集的基函数：1111,1()1,20,iiiiiiiiixxxxxhxxxxxxh在别处1,2,,1in（J.14）11,()0,nnnnnxxxxxhx在别处（J.15）图J.4.线性元的基函数显然，任一hhuV可以表为1()()nhiiiuxux（J.16）有限元方程将变分方程（9）局限在有限元空间上考虑，就得到有限元方程：求有限元解hhuV满足),(),(hhhvfvuahhVv（J.17）注意到hu和hv都可以表示成（J.16）形式，容易看出（J.17）等价于如下的线性方程组：求节点上的近似解1,,nuu满足1(,)(,),1,,nijijiaufjn（J.18）这个线性方程组是三对角的，可以用追赶法求解。可以把微分方程（J.1）、变分方程（J.12）和有限元方程（J.18）比喻为确定“好人”的三种标准：他每时每刻表现都好；大家都说他好；一个遴选委员会说他好。误差估计可以证明，微分方程（J.1）的解u和有限元方程（J.18）的解hu之间的误差满足||||||||||||hhuuuuChu（J.19）其中C是一个常数；||||表示如下定义的2()LI范数：122(,)bavvvvdx（J.20）二维椭圆方程有限元法以二维区域上的Poisson方程第一边值问题为例：2222(,)uufxyxy，(,)xyG（J.21a）|0u（J.21b）其中G是以为边界的一个二维区域。利用Green公式，容易推出相应的变分方程：求10()uHG满足(,)(,)auvfv，10()vHG（J.22）其中函数集合10()HG由满足以下条件的所有函数组成：在边界上为零，且本身及其广义偏导数在区域G上勒贝格可积；(,)dxdyGwvwv（J.23）(,)()Gwvwvawvdxdyxxyy（J.24）二维区域上最常用剖分是形如下图的三角剖分：我们可以相应地构造三角剖分上的线性元。对内点集合hG（例如上图中3，6，5这三个点）中每个节点i，定义其基函数(,)ixy为一个分片线性函数，它在节点i取值为1，而在所有其他节点为0。这样，有限元空间hV中任一元素就可以表示成()()hhiiiGuxux。把它带入到变分方程（J.22）便得有限元方程：求hG上的近似解iu满足(,)(,),hijijhiGaufjG（J.25）高次元可以从两个途径来提高有限元法的精度，一个是加密网格，另一个是利用高次元。例如对于一维问题，可以使用所谓Hermite三次元，它在每一个单元1,iiiIxx上是一个三次多项式，由两个端点上的函数值和导数值总共4个待定参数确定。这时，相应于（J.19）我们有误差估计330||||||||||||khhkkdududuuuChdxdxdx（J.26）对于二维问题也可以使用高次元，但是其定义稍微复杂一点。习题1设边值条件为(0)0u，0)1(u，步长为h=0.25。写出相应的线性元的各个基函数，并图示