微分方程的数值解法课件

第三部分偏微分方程的有限元方法,一边值问题的变分原理,1引论,求使得泛函,达到最大的函数。,（1）等周问题,在长度一定的所有平面封闭曲线中，求所围面积为最大的曲线。,定义：当求泛函在一个函数集合K中的极小（或极大）问题，则该问题称为变分问题。,变分问题与微分方程的定解问题有一定的联系。,（2）初等变分原理,一元二次函数的变分原理,考察J(x)的极值情况。,变分原理：,设,求，使,与求解方程Lx=f等价。,对称正定,多元二次函数的变分原理,求J(x)取极小值的驻点，其中,设,设,则J(x)可表示为：,变分原理：,设矩阵A对称正定，则下列两个命题等价：,(b)是方程的解,上述两个例子表明：,其中,求二次函数的极小值问题和求线性代数方程（组）的解是等价的。,（1）弦平衡的平衡原理与极小位能原理,2两点边值问题的变分原理,考察一根长为l的弦，两端固定在点A(0,0)和B(l,0)。当没有外力作用时，它的位置沿水平方向与X轴重合。设有强度为f(x)的外荷载垂直向下作用在弦上，于是弦发生形变。假定荷载很小，因而发生的形变也很小。用u(x)表示在荷载f(x)的作用下弦的平衡位置。,求弦的平衡位置归结为求解两点边值问题：,设弦处于某一位置u=u(x)，可得到其总位能为,极小位能原理：,其中T是弦的张力。,平衡原理,弦的平衡位置(记为)将在满足边值条件u(0)=0，u(l)=0的一切可能位置中，使位能取极小值。,弦的平衡位置是下列变分问题的解,在数学上，要将某个微分方程的定解问题转化为一个变分问题求解，必须针对已给的定解问题构造一个相应的泛函，并证明定解问题的解与泛函极值问题的解等价。,有限元方法正是利用这种等价性（边值问题与变分问题的等价性），先将微分方程定解问题转化为变分问题（或变分方程）的求解问题，然后再设法近似求解变分问题（或变分方程）。,（2）两点边值问题的变分原理,构造泛函,考察二阶常微分方程边值问题：,引入泛函算子,则,变分问题,与前述二阶常微分方程边值问题相应的变分问题是,其中,求，使,变分原理（变分问题与边值问题的等价性）,设，是边值问题,的解，则使J(u)达到极小值；,反之，若使J(u)达到极小值，则是边值问题的解。,其中,是强制边界条件，是自然边界条件，区别这两类边界条件在用有限元方法求解边值问题时很重要。,（3）虚功原理,对两点边值问题：,其中,虚功原理,，且满足变分方程：,设，以v乘方程两端，沿a,b积分，并利用，得变分方程,对任意,在力学里，表示虚功,设，则是边值问题解的充要条件是：,对于复杂的边界条件，边值问题的求解一般是困难的。若将微分方程化为相应的变分问题或变分方程，则只需处理强加边界条件，无需处理自然边界条件（自然边界条件已包含于变分问题中泛函的构造或已包含于给出的变分方程之中）。这一特点对研究微分方程离散化方法及其数值解带来了极大的方便。,3二阶椭圆边值问题的变分原理,（1）极小位能原理,模型方程,其中G是平面有界区域。,构造泛函,引入泛函算子,则,变分问题,与前述二阶椭圆边值问题相应的变分问题是,求，使,其中,变分原理（变分问题与边值问题的等价性）,对第一边值问题，无论齐次或非齐次边界条件，泛函是一样的，只是边界条件要作为强加边值条件加在所取的函数类上。,设，是二阶椭圆边值问题的解，则使J(u)达到极小值；,反之，若使J(u)达到极小值，则是二阶椭圆边值问题的解。,其中,对第二、三类边值问题，无论齐次或非齐次边界条件，二次泛函形式相对于第一边值问题有所改变，但函数类的选取与边界条件无关。,（2）虚功原理,问题,其中,设，以v乘方程两端后在G上积分，并利用Green公式，得变分方程,虚功原理,在力学里，表示虚功,设是边值问题的解，则对任意，满足变分方程。,反之，若，且对任意满足变分方程，则为边值问题的解。,与极小位能原理类似，第一类边界条件为强加边界条件，第二、三类边界条件为自然边界条件。,虚功原理比极小位能原理应用更广。,目的：求解相应的变分问题或相应的变分方程。,Ritz方法是近似求解变分问题(即二次泛函极小值)的算法。Galerkin方法是近似求解变分方程的算法，这两种算法统称为Ritz-Galerkin方法。,Ritz-Galerkin方法的基本思想,以下用V表示等Sobolev空间，L表示微分算子，(u,v)为由L及边值条件决定的双线性泛函。,4Ritz-Galerkin方法,用有限维空间的函数代替变分问题(或变分方程)中无限维空间的函数，从而在有限维函数空间中求变分问题(或变分方程)的近似解，并要求当有限维空间的维数不断增加时，有限维近似解逼近原变分问题(或变分方程)的解。,由极小位能原理得出的变分问题为：,Ritz方法：求变分问题的近似解。,（1）Ritz方法,求，使,其中，,设是V的n维子空间，是的一组基底(称为基函数)。中任一元素可表示为,即选择适当的，使取极小值。,求，使,Ritz方法：,展开,令,则满足,解出代入，则得,Ritz方法步骤为：,根据最小位能原理构造相应于微分方程或物理问题的变分问题；,取作为的一组基底，即用近似代替无穷维空间V；,求解关于的线性代数方程组。,由虚功原理得出的变分方程为：,Galerkin方法：求变分方程的近似解。,（2）Galerkin方法,设是V的n维子空间，是的一组基底(称为基函数)。中任一元素可表示为,即选择适当的，使取极小值。,Galerkin方法：,求，使对，满足,由的任意性，取作为v，则得,将代入变分方程，则,解出代入，则得,Galerkin步骤为：,根据虚功原理构造相应于微分方程或物理问题的变分方程；,取作为的一组基底，即用近似代替无穷维空间V；,求解关于的线性代数方程组。,取作为v，将代入变分方程，得到满足的方程组：,有限元法广泛应用的原因,Ritz-Galerkin方法应用的困难,基函数选取必须满足强加边界条件，因此选取困难；,计算量、存储量巨大；,方程组求解病态严重。,充分发挥了变分形式和Ritz-Galerkin方法的优点；,摆脱了传统的基函数取法；,各种问题的结构程序格式统一。,有限元方法基于变分原理，又具有差分方法的一些特点，并且适于较复杂的区域和不同粗细的网格。,二椭圆型方程的有限元方法,差分法解偏微分方程，解得的结果就是准确解u在节点上的近似值；,Ritz-Galerkin方法得到近似的解析解，但对一般区域，却往往难以实现。,有限元方法与传统Ritz-Galerkin方法的差别在于有限维函数空间的构造方法。Ritz-Galerkin方法选用的基函数在整个定解区域上整体光滑，有限元则取分段或分片连续且局部非零的基函数。,考虑两点边值问题：,1一维问题的线性元,将区间a,b分割为n个子区间。,第i个单元记为，其长度。,（1）试探函数与试探函数空间,设,则称为试探函数空间，称为试探函数。,（2）用单元形状函数表示试探函数,设在节点上试探函数在节点上的一组值为,最简单的试探函数空间由分段线性函数组成。,在第i个单元上的线性插值函数为,即,当时，的（线性）插值公式称为（线性）单元形状函数。,把每个单元形状函数合并起来，就得到整个区间a,b上都有定义的函数：,为使分段插值标准化，通常用仿射变换,显然,把变到，令,则,变为,或,定义基函数系,（3）用节点基函数表示试探函数,线性无关，它们可组成试探函数空间的基，常称为节点基函数。,几何形状如图,任一试探函数可表示为,用这类插值型基函数，可以构造出适合各种边界条件的试探函数。,若借助前述放射变换,节点基函数可用变量表示为,直接形成有限元方程,(a)把表达式代入泛函；,（4）从Ritz方法出发形成有限元方程,(b)将泛函表达式中积分区间a,b变到0,1；,(c)由达到极小值的条件,得到含的有限元方程,这儿,(d)解出有限元方程的数值解，就得到使二次泛函取极小的近似函数(有限元解),有限元方程可用矩阵表示为,其中,称为总刚矩阵。,工程中形成有限元方程时，通常先在每个单元上形成单元矩阵（称为单元刚度矩阵），然后由单元刚度矩阵形成总刚度矩阵（称为总体合成）。,用单元刚度分析形成有限元方程,(a)把按单元组织，则在第i个单元上，令,其中称为单元刚度矩阵。各元素可计算得到。,再把扩展成nn矩阵，使其第i1行、第i行和第i1列、第i列交叉位置的元素就是单元刚度矩阵的四个元素，其余全为零（只是第一行，第一列元素非零）。即,记,则,其中称为总刚矩阵。,(b)由达到极小值的条件,(c)解出有限元方程的数值解，就得到使二次泛函取极小的近似函数(有限元解),得到有限元方程。,（5）从Galerkin方法出发形成有限元方程,把表达式代入变分方程,对前面的两点边值问题，变分方程变为,其中,与Ritz方法相比，Galerkin方法形成的有限元方程其系数矩阵就是总刚矩阵。,该方程即为Galerkin法形成的有限元方程。,由Galerkin方法推导有限元方程更加方便直接，且适用面广。,若希望在每个单元上提高逼近的精确度，则可通过提高插值多项式次数来实现，,在单元上可构造一、二、三及高次插值多项式，其方法有两种：,2一维问题的高次元,整个问题计算的全过程除分析单元插值外，均与前面框架类似。,Lagrange型：在单元内部增加一些插值节点。,Hermite型：在节点引进一阶、二阶乃至更高阶导数。,线性元(Lagrange型),要求：在每一个单元上是一次多项式，在单元节点处连续。插值条件：在单元的两个端点取指定值。,二次元(Lagrange型),要求：在每一个单元上是二次多项式，在单元节点处连续。插值条件：在单元的两个端点及单元中点取指定值。,三次元（Hermite型）,要求：在每一个单元上是三次多项式，在单元节点处连续。插值条件：在两个端点取指定的函数值和一阶导数值。,采用高次元，有限元方程形成的方法和线性元类似，但工作量增加。一是计算积分的复杂性增加，二是矩阵的带宽增加。,高次元的主要优点是收敛阶高，且提高了函数逼近的光滑性。,假定区域G可以分割成有限个矩形的和，且每个小矩形（单元）的边和坐标轴平行。,3二维问题的矩形元,通过仿射变换,采用矩形剖分后，任一个矩形,总可变成单位正方形,如果在上造出单元形状函数，就可得到试探函数。而上的形状函数可通过先在上造出形状函数，再通过仿射变化而得到。,在上构造形状函数，也采用Lagrange型和Hermite型插值。,Lagrange型：根据若干插值节点处的函数值决定插值函数。,Hermite型：根据若干插值节点处的函数值、一阶偏导数乃至更高阶偏导数决定插值函数。,（1）Lagrange型公式,双一次插值,插值条件：给定顶点上的函数值,求：双线性函数,满足,设,令,由为双线性函数，可求得,令,则,通过仿射变换消去、，就得到上的形状函数。把这些函数按单元叠加，即对所有单元求和，就得到G上的试探函数。,实际计算时，并不消去中间变量、，因为计算刚度矩阵元素（定积分）用、作自变量更为方便。,插值条件：给定II上九个插值节点(0,0)、(1/2,0)、(1,0)、(0,1/2)、(1/2,1/2)、(1,1/2)、(0,1)、(1/2,1)、(1,1)的函数值。,求：双二次函数,满足,双二次插值,故,通过仿射变换消去、，就得到上的形状函数。,令,由为二次函数，可求得,设,插值条件：给定II上十六个插值节点(见图)。,求：双三次函数,满足,设,双三次插值,故,令,由为三次函数，可求得,可以在四个顶点分别给定函数值、两个一阶偏导数的值和二阶混合偏导数的值(共十六条件)，确定一个双三次多项式的十六个系数。,（2）Hermite型公式,Lagrange型公式中不出现导数，这样的试探函数只属于。为了得到属于的试探函数，需要Hermite型插值公式。,双三次多项式含有十六项：,简单且常用的是不完全的双三次多项式插值。它去掉双三次多项式中的项。,插值条件：给定II上四个插值节点。,求：不完全双三次函数,满足,四个顶点处的函数值等于在该点的函数值；四个顶点处的值等于在该点的值；四个顶点处的值等于在该点的值。,根据仿射变换,则可将原插值问题转化为II上的插值问题。,满足,四个顶点处的函数值等于在该点的函数值；四个顶点处的值等于在该点的值乘以x；四个顶点处的值等于在该点的值乘以y。,插值条件：给定II上四个插值节点(0,0)、(1,0)、(0,1)、(1,1)。,求：不完全双三次函数,类似于Lagrange型公式的构造，可以求得上的形状函数。,在三角形元的有限元方法中，先将定解区域G化分为若干个小三角形（称作单元）。然后在每个单元上构造插值型函数，并用分片函数（但整体连续的函数）代替变分问题或变分方程中所需求解的函数。,4二维问题的三角形元,用有限元求解二维椭圆边值问题时，应用最广的是三角形元。,（1）三角剖分,将定解区域化分成若干个小三角形单元时应注意：,为了保证有限元解的精确度和收敛性，并避免其离散后代数方程组系数矩阵的病态性，网格剖分中疏密的过渡不要太陡。,错误,为了保证有限元解有较好的精度，每个单元中应尽量避免出现大的钝角。,应避免,单元顶点的编号顺序可以任意，但节点编号顺序将影响有限元方程组系数矩阵的结构（带宽）。,为了方便构造插值型函数，要求每个单元的顶点是相邻单元的顶点。,（2）面积坐标及有关公式,在三角形单元上构造插值型函数，并不简单类同于矩形单元。,面积坐标,考虑一个面积为S的三角形单元，其顶点按反时针顺序记为i,j,k。在此单元内部任取一点p(x,y)，连接p和三个顶点，此单元则被分成三个小三角形它们的面积记为和。,记,单元内任一点p(x,y)的位置与三维数一一对应，称为面积坐标。,面积坐标和直角坐标之间的关系：,面积坐标与直角坐标的关系,面积坐标与坐标系无关，这是采用面积坐标的优点,面积坐标有性质：,由于面积坐标满足，将其代入,得：,任意三角形到标准等腰直角三角形的变换,将看作是某一平面的坐标，则上式表明了一种变换。可以证明它把X-Y坐标系的任意三角形单元映射为坐标系中的标准等腰直角三角形单元。且,即该变换是仿射变换。,这个变换的Jacobi行列式,该变换除了能将三角单元仿射变换为标准三角单元，还能将三角单元上的插值型函数变换为标准三角单元上的同类型函数。因此，采用面积坐标可使计算工作简单化、标准化。,另外，利用面积坐标表示的齐次多项式在(i,j,k)上的积分也非常容易计算。即,其中p,q,r是任意非负整数。,一般在三角形元上，构造一个m次完全多项式,（3）Lagrange型公式,两个变量x,y的高次多项式可用Pascal三角形表示：,故Lagrange型插值需要个节点的函数插值条件。,逼近u(x,y)时，该多项式具有的项数为,一次多