微分方程的数值解法

微分方程的数值解法第一页，共六十二页，编辑于2023年，星期六定义：当求泛函在一个函数集合K中的极小（或极大）问题，则该问题称为变分问题。

变分问题与微分方程的定解问题有一定的联系。（2）初等变分原理①一元二次函数的变分原理考察J(x)的极值情况。变分原理：设求，使与求解方程Lx=f

等价。第二页，共六十二页，编辑于2023年，星期六对称正定②多元二次函数的变分原理求J(x)取极小值的驻点，其中设设则J(x)可表示为：第三页，共六十二页，编辑于2023年，星期六变分原理：设矩阵A对称正定，则下列两个命题等价：求，使(a)(b)

是方程的解上述两个例子表明：其中

求二次函数的极小值问题和求线性代数方程（组）的解是等价的。第四页，共六十二页，编辑于2023年，星期六（1）弦平衡的平衡原理与极小位能原理2两点边值问题的变分原理

考察一根长为l的弦，两端固定在点A(0,0)和B(l,0)。当没有外力作用时，它的位置沿水平方向与X轴重合。设有强度为f(x)的外荷载垂直向下作用在弦上，于是弦发生形变。假定荷载很小，因而发生的形变也很小。用u(x)表示在荷载f(x)的作用下弦的平衡位置。第五页，共六十二页，编辑于2023年，星期六求弦的平衡位置归结为求解两点边值问题：设弦处于某一位置u=u(x)，可得到其总位能为

极小位能原理：其中T是弦的张力。平衡原理

弦的平衡位置(记为)将在满足边值条件u(0)=0，u(l)=0的一切可能位置中，使位能取极小值。弦的平衡位置是下列变分问题的解第六页，共六十二页，编辑于2023年，星期六

在数学上，要将某个微分方程的定解问题转化为一个变分问题求解，必须针对已给的定解问题构造一个相应的泛函，并证明定解问题的解与泛函极值问题的解等价。

有限元方法正是利用这种等价性（边值问题与变分问题的等价性），先将微分方程定解问题转化为变分问题（或变分方程）的求解问题，然后再设法近似求解变分问题（或变分方程）。第七页，共六十二页，编辑于2023年，星期六（2）两点边值问题的变分原理①构造泛函考察二阶常微分方程边值问题：引入泛函算子则第八页，共六十二页，编辑于2023年，星期六②变分问题

与前述二阶常微分方程边值问题相应的变分问题是其中求，使第九页，共六十二页，编辑于2023年，星期六③变分原理（变分问题与边值问题的等价性）设，是边值问题的解，则使J(u)

达到极小值；

反之，若使J(u)

达到极小值，则是边值问题的解。其中

是强制边界条件，是自然边界条件，区别这两类边界条件在用有限元方法求解边值问题时很重要。第十页，共六十二页，编辑于2023年，星期六（3）虚功原理对两点边值问题：其中虚功原理，且满足变分方程：

设，以v乘方程两端，沿[a,b]积分，并利用，得变分方程对任意在力学里，表示虚功设，则

是边值问题解的充要条件是：第十一页，共六十二页，编辑于2023年，星期六

对于复杂的边界条件，边值问题的求解一般是困难的。若将微分方程化为相应的变分问题或变分方程，则只需处理强加边界条件，无需处理自然边界条件（自然边界条件已包含于变分问题中泛函的构造或已包含于给出的变分方程之中）。这一特点对研究微分方程离散化方法及其数值解带来了极大的方便。第十二页，共六十二页，编辑于2023年，星期六3二阶椭圆边值问题的变分原理（1）极小位能原理模型方程其中G是平面有界区域。①构造泛函引入泛函算子则第十三页，共六十二页，编辑于2023年，星期六②变分问题与前述二阶椭圆边值问题相应的变分问题是求，使其中第十四页，共六十二页，编辑于2023年，星期六③变分原理（变分问题与边值问题的等价性）

对第一边值问题，无论齐次或非齐次边界条件，泛函是一样的，只是边界条件要作为强加边值条件加在所取的函数类上。

设，是二阶椭圆边值问题的解，则使J(u)达到极小值；

反之，若使J(u)

达到极小值，则是二阶椭圆边值问题的解。其中

对第二、三类边值问题，无论齐次或非齐次边界条件，二次泛函形式相对于第一边值问题有所改变，但函数类的选取与边界条件无关。第十五页，共六十二页，编辑于2023年，星期六（2）虚功原理问题其中

设，以v乘方程两端后在G上积分，并利用Green公式，得变分方程第十六页，共六十二页，编辑于2023年，星期六虚功原理在力学里，表示虚功

设是边值问题的解，则对任意，满足变分方程。

反之，若，且对任意满足变分方程，则为边值问题的解。

与极小位能原理类似，第一类边界条件为强加边界条件，第二、三类边界条件为自然边界条件。

虚功原理比极小位能原理应用更广。第十七页，共六十二页，编辑于2023年，星期六目的：求解相应的变分问题或相应的变分方程。

Ritz方法是近似求解变分问题(即二次泛函极小值)的算法。Galerkin方法是近似求解变分方程的算法，这两种算法统称为Ritz-Galerkin方法。Ritz-Galerkin方法的基本思想

以下用V表示等Sobolev空间，L表示微分算子，(u,v)为由L及边值条件决定的双线性泛函。4Ritz-Galerkin方法

用有限维空间的函数代替变分问题(或变分方程)中无限维空间的函数，从而在有限维函数空间中求变分问题(或变分方程)的近似解，并要求当有限维空间的维数不断增加时，有限维近似解逼近原变分问题(或变分方程)的解。第十八页，共六十二页，编辑于2023年，星期六由极小位能原理得出的变分问题为：Ritz方法：求变分问题的近似解。（1）Ritz方法求，使其中，

设是V

的n维子空间，是的一组基底(称为基函数)。中任一元素可表示为即选择适当的，使取极小值。求，使Ritz方法：第十九页，共六十二页，编辑于2023年，星期六展开令则满足解出代入，则得第二十页，共六十二页，编辑于2023年，星期六Ritz方法步骤为：

根据最小位能原理构造相应于微分方程或物理问题的变分问题；

取作为的一组基底，即用近似代替无穷维空间V；

根据二次函数取极值的必要条件，得到中所满足的方程组：

求解关于的线性代数方程组。第二十一页，共六十二页，编辑于2023年，星期六由虚功原理得出的变分方程为：Galerkin方法：求变分方程的近似解。（2）Galerkin方法

设是V

的n维子空间，是的一组基底(称为基函数)。中任一元素可表示为即选择适当的，使取极小值。Galerkin方法：求，使对，满足第二十二页，共六十二页，编辑于2023年，星期六由的任意性，取作为v

，则得将代入变分方程，则解出代入，则得第二十三页，共六十二页，编辑于2023年，星期六Galerkin步骤为：

根据虚功原理构造相应于微分方程或物理问题的变分方程；

取作为的一组基底，即用近似代替无穷维空间V；

求解关于的线性代数方程组。取作为v

，将代入变分方程，得到满足的方程组：第二十四页，共六十二页，编辑于2023年，星期六有限元法广泛应用的原因Ritz-Galerkin方法应用的困难①基函数选取必须满足强加边界条件，因此选取困难；②计算量、存储量巨大；③方程组求解病态严重。

充分发挥了变分形式和Ritz-Galerkin方法的优点；②摆脱了传统的基函数取法；③各种问题的结构程序格式统一。第二十五页，共六十二页，编辑于2023年，星期六

有限元方法基于变分原理，又具有差分方法的一些特点，并且适于较复杂的区域和不同粗细的网格。二椭圆型方程的有限元方法

差分法解偏微分方程，解得的结果就是准确解u在节点上的近似值；Ritz-Galerkin方法得到近似的解析解，但对一般区域，却往往难以实现。

有限元方法与传统Ritz-Galerkin方法的差别在于有限维函数空间的构造方法。Ritz-Galerkin方法选用的基函数在整个定解区域上整体光滑，有限元则取分段或分片连续且局部非零的基函数。第二十六页，共六十二页，编辑于2023年，星期六考虑两点边值问题：1一维问题的线性元

将区间[a,b]分割为n个子区间。第i个单元记为，其长度。（1）试探函数与试探函数空间设则称为试探函数空间，称为试探函数。第二十七页，共六十二页，编辑于2023年，星期六（2）用单元形状函数表示试探函数设在节点上试探函数在节点上的一组值为最简单的试探函数空间由分段线性函数组成。在第i个单元上的线性插值函数为即

当时，的（线性）插值公式称为（线性）单元形状函数。第二十八页，共六十二页，编辑于2023年，星期六

把每个单元形状函数合并起来，就得到整个区间[a,b]上都有定义的函数：第二十九页，共六十二页，编辑于2023年，星期六为使分段插值标准化，通常用仿射变换显然把变到，令则变为或第三十页，共六十二页，编辑于2023年，星期六定义基函数系（3）用节点基函数表示试探函数第三十一页，共六十二页，编辑于2023年，星期六

线性无关，它们可组成试探函数空间的基，常称为节点基函数。几何形状如图ab任一试探函数可表示为

用这类插值型基函数，可以构造出适合各种边界条件的试探函数。第三十二页，共六十二页，编辑于2023年，星期六若借助前述放射变换节点基函数可用变量表示为第三十三页，共六十二页，编辑于2023年，星期六①直接形成有限元方程(a)把表达式代入泛函；（4）从Ritz方法出发形成有限元方程(b)将泛函表达式中积分区间[a,b]变到[0,1]；(c)由达到极小值的条件得到含的有限元方程这儿(d)解出有限元方程的数值解，就得到使二次泛函取极小的近似函数(有限元解)第三十四页，共六十二页，编辑于2023年，星期六有限元方程可用矩阵表示为其中称为总刚矩阵。第三十五页，共六十二页，编辑于2023年，星期六

工程中形成有限元方程时，通常先在每个单元上形成单元矩阵（称为单元刚度矩阵），然后由单元刚度矩阵形成总刚度矩阵（称为总体合成）。②用单元刚度分析形成有限元方程(a)把按单元组织，则在第i个单元上，令其中称为单元刚度矩阵。各元素可计算得到。第三十六页，共六十二页，编辑于2023年，星期六

再把扩展成nn矩阵，使其第i1行、第i行和第i1列、第i列交叉位置的元素就是单元刚度矩阵的四个元素，其余全为零（只是第一行，第一列元素非零）。即记则其中称为总刚矩阵。第三十七页，共六十二页，编辑于2023年，星期六(b)由达到极小值的条件(c)解出有限元方程的数值解，就得到使二次泛函取极小的近似函数(有限元解)得到有限元方程。第三十八页，共六十二页，编辑于2023年，星期六（5）从Galerkin方法出发形成有限元方程把表达式代入变分方程对前面的两点边值问题，变分方程变为其中

与Ritz方法相比，Galerkin方法形成的有限元方程其系数矩阵就是总刚矩阵。该方程即为Galerkin法形成的有限元方程。

由Galerkin方法推导有限元方程更加方便直接，且适用面广。第三十九页，共六十二页，编辑于2023年，星期六

若希望在每个单元上提高逼近的精确度，则可通过提高插值多项式次数来实现，

在单元上可构造一、二、三及高次插值多项式，其方法有两种：2一维问题的高次元

整个问题计算的全过程除分析单元插值外，均与前面框架类似。①Lagrange型：在单元内部增加一些插值节点。②Hermite型：在节点引进一阶、二阶乃至更高阶导数。第四十页，共六十二页，编辑于2023年，星期六①线性元(Lagrange型)要求：在每一个单元上是一次多项式，在单元节点处连续。插值条件：在单元的两个端点取指定值。②二次元(Lagrange型)要求：在每一个单元上是二次多项式，在单元节点处连续。插值条件：在单元的两个端点及单元中点取指定值。③三次元（Hermite型）要求：在每一个单元上是三次多项式，在单元节点处连续。插值条件：在两个端点取指定的函数值和一阶导数值。第四十一页，共六十二页，编辑于2023年，星期六

采用高次元，有限元方程形成的方法和线性元类似，但工作量增加。一是计算积分的复杂性增加，二是矩阵的带宽增加。

高次元的主要优点是收敛阶高，且提高了函数逼近的光滑性。第四十二页，共六十二页，编辑于2023年，星期六

假定区域G可以分割成有限个矩形的和，且每个小矩形（单元）的边和坐标轴平行。3二维问题的矩形元通过仿射变换采用矩形剖分后，任一个矩形总可变成单位正方形

如果在上造出单元形状函数，就可得到试探函数。而上的形状函数可通过先在上造出形状函数，再通过仿射变化而得到。第四十三页，共六十二页，编辑于2023年，星期六

在上构造形状函数，也采用Lagrange型和Hermite型插值。Lagrange型：根据若干插值节点处的函数值决定插值函数。Hermite型：根据若干插值节点处的函数值、一阶偏导数乃至更高阶偏导数决定插值函数。第四十四页，共六十二页，编辑于2023年，星期六（1）Lagrange型公式①双一次插值插值条件：给定顶点上的函数值求：双线性函数满足设令由为双线性函数，可求得第四十五页，共六十二页，编辑于2023年，星期六令则

通过仿射变换消去、

，就得到上的形状函数。把这些函数按单元叠加，即对所有单元求和，就得到G上的试探函数。

实际计算时，并不消去中间变量、，因为计算刚度矩阵元素（定积分）用、作自变量更为方便。第四十六页，共六十二页，编辑于2023年，星期六插值条件：给定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)的函数值。求：双二次函数满足②双二次插值第四十七页，共六十二页，编辑于2023年，星期六故

通过仿射变换消去、

，就得到上的形状函数。令由为二次函数，可求得设第四十八页，共六十二页，编辑于2023年，星期六插值条件：给定II上十六个插值节点(见图)。求：双三次函数满足设③双三次插值第四十九页，共六十二页，编辑于2023年，星期六故令由为三次函数，可求得第五十页，共六十二页，编辑于2023年，星期六

可以在四个顶点分别给定函数值、两个一阶偏导数的值和二阶混合偏导数的值(共十六条件)，确定一个双三次多项式的十六个系数。（2）Hermite型公式Lagrange型公式中不出现导数，这样的试探函数只属于。为了得到属于的试探函数，需要Hermite型插值公式。双三次多项式含有十六项：

简单且常用的是不完全的双三次多项式插值。它去掉双三次多项式中的项。第五十一页，共六十二页，编辑于2023年，星期六插值条件：给定II上四个插值节点。求：不完全双三次函数满足四个顶点处的函数值等于在该点的函数值；四个顶点处的值等于在该点的值；四个顶点处的值等于在该点的值。根据仿射变换则可将原插值问题转化为II上的插值问题。第五十二页，共六十二页，编辑于2023年，星期六满足四个顶点处的函数值等于在该点的函数值；四个顶点处的值等于在该点的值乘以x；四个顶点处的值等于在该点的值乘以y。插值条件：给定II上四个插值节点(0,0)、(1,0)、(0,1)、(1,1)。求：不完全双三次函数

类似于Lagrange型公式的构造，可以求得上的形状函数。第五十三页，共六十二页，编辑于2023年，星期六

在三角形元的有限元方法中，先将定解区域G化分为若干个小三角形（称作单元）。然后在每个单元上构造插值型函数，并用分片函数（但整体连续的函数）代替变分问题或变分方程中所需求解的函数。4二维问题的三角形元

用有限元求解二维椭圆边值问题时，应用最广的是三角形元。第五十四页，共六十二页，编辑于2023年，星期六（1）三角剖分

将定解区域化分成若干个小三角形单元时应注意：③为了保证有限元解的精确度和收敛性，并避免其离散后代数方程组系数矩阵的病态性，网格剖分中疏密的过渡不要太陡。错误为了保证有限元解有较好的精度，每个单元中应尽量避免出现大的钝角。应避免④单元顶点的编号顺序可以任意，但节点编号顺序将影响有限元方程组系数矩阵的结构（带宽）。①为了方便构造插值型函数，要求每个单元的顶点是相邻单元的顶点。第五十五页，共六十二页，编辑于2023年，星期六（2）面积坐标及有关公式

在三角形单元上构造插值型函数，并不简单类同于矩形单元。①面积坐标

考虑一个面积为S的三角形单元，其顶点按反时针顺序记为i,j,k。在此单元内部任取一点p(x,y)，连接p和三个顶点，此单元则被分成三个小三角形它们的面积记为和。ijkP(x,y)记单元内任一点p(x,y)的位置与三维数一一对应，称