今天介绍关于计算机求解微分方程的六大数值计算的方法。
1.有限元法
有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式,借助于变分原理或加权余量法,将微分方程离散求解。采用不同的权函数和插值函数 形式,便构成不同的有限元方法。
在有限元方法中,把计算域离散剖分为有限个互不重叠且相互连接的单元,在每个单元内选择基函数,用单元基函数的线形组合来逼近单元中的真解,整个计算域上总体的基函数可以看为由每个单元基函数组成的,则整个计算域内的解可以看作是由所有单元 上的近似解构成。
根据所采用的权函数和插值函数的不同 ,有限元方法也分为多种计算格式。从权函数的选择来说,有配置法、矩量法、最小二乘法和伽辽金法,从计算单元网格的形状来划分,有三角形网格、四边形网格和多边形网格,从插值函数的精度来划分,又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。
对于权函数,伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ;最小二乘法是令权函数等于余量本身,而内积的极小值则为对代求系数的平方误差最小;在配置法中,先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程,即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成,但也有采用三角函数或指数函数组成的乘积表示,但最常用的多项式插值函数。有限元插值函数分为两大类,一类只要求插值多项式本身在插值点取已知值,称为拉格朗日(Lagrange)多项式插值;另一种不仅要求插值多项式本身,还要求它的导数值在插值点取已知值,称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标,有对称和不对称等。常采用的无因次坐标是一种局部坐标系,它的定义取决于单元的几何形状,一维看作长度比,二维看作面积比,三维看作体积比。在二维有限元中,三角形单元应用的最早,近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元,常采用的插值函数为有La grange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线 性插值函数、二阶或更高阶插值函数等。
对于有限元方法,其基本思路和解题步骤可归纳为
- 建立积分方程,根据变分原理或方程余量与权函数正交化原理,建立与微分方程初边值问题等价的积分表达式,这是有限元法的出发点。
- 区域单元剖分,根据求解区域的形状及实际问题的物理特点,将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作,这部分工作量比较大,除了给计算单元和节点进行编号和确定相互之间的关系之外,还要表示节点的位置坐标,同时还需要列出自然边界和本质边界的节点序号和相应的边界值。
- 确定单元基函数,根据单元中节点数目及对近似解精度的要求,选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的,由于各单元 具有规则的几何形状,在选取基函数时可遵循一定的法则。
- 单元分析:将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近;再将近似函数代入积分方程,并对单元区域进行积分,可获得含有待定系数(即单元中各节点 的参数值)的代数方程组,称为单元有限元方程。
- 总体合成:在得出单元有限元方程之后,将区域中所有单元有限元方程按一定法则进 行累加,形成总体有限元方程。
- 边界条件的处理:一般边界条件有三种形式,分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件,一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件,需按一定法 则对总体有限元方程进行修正满足。
- 解有限元方程:根据边界条件修正的总体有限元方程组,是含所有待定未知量的封闭方程组,采用适当的数值计算方法求解,可求得各节点的函数值。
2. 多重网格方法
多重网格方法通过在疏密不同的网格层上进行迭代,以平滑不同频率的误差分量.具有收敛速度快,精度高等优点.
多重网格法基本原理微分方程的误差分量可以分为两大类,一类是频率变化较缓慢的低频分量;另一类是频率高,摆动快的高频分量。一般的迭代方法可以迅速地将摆动误差衰减,但对那些低频分量,迭代法的效果不是很显著。高频分量和低频分量是相对的,与网格尺度有关,在细网格上被视为低频的分量,在粗网格上可能为高频分量。
多重网格方法作为一种快速计算方法,迭代求解由偏微分方程组离散以后组成的代数方程组,其基本原理在于一定的网格最容易消除波长与网格步长相对应的误差分量。该方法采用不同尺度的网格,不同疏密的网格消除不同波长的误差分量,首先在细网格上采用迭代法,当收敛速度变缓慢时暗示误差已经光滑,则转移到较粗的网格上消除与该层网格上相对应的较易消除的那些误差分量,这样逐层进行下去直到消除各种误差分量,再逐层返回到细网格上。目前两层网格方法从理论上已证明是收敛的,并且其收敛速度与网格尺度无关。
多重网格法是迭代法与粗网格修正的组合,经过证明迭代法可迅速地将那些高频分量去掉,粗网格修正则可以帮助消除那些光滑了的低频分量,而对那些高频分量基本不起作用。科研中国SciEi.com 在多重网格计算中,需要一些媒介把细网格上的信息传递到粗网格上去,同时还需要一些媒介把粗网格上的信息传递到细网格上去。限制算子Iih(i-1)h是把细网格i-1层上的残余限制到粗网格i层上的算子,最简单的算子是平凡单射,另外还有特殊加权限制;插值算子Iih(i-1)h是把粗网格i层上的结果插值到细网格i-1层上的算子,一般采用线性插值或完全加权限制算子。
3.有限差分方法
有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑,可分为中心格式和逆风格式。考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分方法主要适用于有结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等,其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。
4.有限体积法
有限体积法(Finite Volume Method)又称为控制体积法。其基本思路是:将计算区域划分为一系列不重复的控制体积,并使每个网格点周围有一个控制体积;将待解的微分方程对每一个控制体积积分,便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分,必须假定值在网格点之间的变化规律,即假设值的分段的分布的分布剖面。从积分区域的选取方法看来,有限体积法属于加权剩余法中的子区域法;从未知解的近似方法看来,有限体积法属于采用局部近似的离散方法。简言之,子区域法属于有限体积发的基本方法。
有限体积法的基本思路易于理解,并能得出直接的物理解释。离散方程的物理意义,就是因变量在有限大小的控制体积中的守恒原理,如同微分方程表示因变量在无限小的控 制体积中的守恒原理一样。限体积法得出的离散方程,要求因变量的积分守恒对任意一组控制体积都得到满足,对整个计算区域,自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法,例如有限差分法,仅当网格极其细密时,离散方程才满足积分守恒;而有限体积法即使在粗网格情况下,也显示出准确的积分守恒。就离散方法而言,有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律(既插值函数),并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值 ,这与有限差分法相类似;但有限体积法在寻求控制体积的积分时,必须假定值在网格点之间的分布,这又与有限单元法相类似。在有限体积法中,插值函数只用于计算控制体积的积分,得出离散方程之后,便可忘掉插值函数;如果需要的话,可以对微分方程 中不同的项采取不同的插值函数。
5.近似求解的误差估计方法
近似求解的误差估计方法共有三大类:单元余量法,通量投射法及外推法.
单元余量法广泛地用于以FEM离散的误差估计之中,它主要是估计精确算子的余量,而不是整套控制方程的全局误差.这样就必须假定周围的单元误差并不相互耦合,误差计算采用逐节点算法进行.单元余量法的各种不同做法主要来自对单元误差方程的边界条件的不同处理办法.基于此,该方法能够有效处理局部的残余量,并能成功地用于网格优化程序.
通量投射法的基本原理来自一个很简单的事实:精确求解偏微分方程不可能有不连续的微分,而近似求解却可以存在微分的不连续,这样产生的误差即来自微分本身,即误差为系统的光滑求解与不光滑求解之差.该方法与单元余量法一样,对节点误差采用能量范数,故也能成功地用于网格优化程序.单元余量法及通量投射法都局限于局部的误差计算(采用能量范数),误差方程的全局特性没有考虑.另外计算的可行性(指误差估计方程的计算时间应小于近似求解计算时间)不能在这两种方法中体现,因为获得的误差方程数量,阶数与流场控制方程相同.
外推法是指采用后向数值误差估计思想由精确解推出近似解的误差值.各类文献中较多地采用Richardson外推方法来估计截断误差.无论是低阶还是高阶格式,随着网格的加密数值计算结果都会趋近于准确解.但由于计算机内存与计算时间的限制,实际上不能采用这种网格无限加密的办法.由Richardson所发展起来的外推方法,可以利用在不同疏密网格上得出的结果估计相应的收敛解,可以估计所用离散方法截断误差的阶数,可以估计所得数值计算的截断误差.该方法有很大的局限性,不能简单地用于复杂湍流流动;并且在数值计算中数值解必须单调地趋近于其收敛值.而文献提出的单网格后向误差估计思想,在采用有限元法FEM,有限容积法FVM时均有应用,并且还用于网格优化程序,但该方法也不能用于复杂湍流流动的数值分析.
6.多尺度计算方法
近年来发展的多尺度计算方法包括均匀化方法、非均匀化多尺度方法、以及小波数值均匀化方法、多尺度有限体积法、多尺度有限元法等。
均匀化方法是一种多尺度分析的方法。该方法通过对单胞问题的求解,把细观尺度的信息映射到宏观尺度上,从而推导出宏观尺度上的均匀化等式,即可在宏观尺度上求解原问题。均匀化方法在很多科学和工程应用中取得了巨大成功,但这种方法建立在系数细观结构周期性假设的基础上,因此应用范围受到了很大限制。
鄂维南等提出的非均匀化多尺度方法,是构造多尺度计算方法的一般框架。该方法有两个重要的组成部分:基于宏观变量的整体宏观格式和由微观模型来估计缺少的宏观数据,多尺度问题的解通过这两部分共同得到。
小波数值均匀化方法是由Dorbonuat、Enguqist提出的求解椭圆型方程的新型方法。该方法基于多分辨分析,在细尺度上建立原方程的离散算子,然后对离散算子进行小波变换,得到了大尺度上的数值均匀化算子。此方法在大尺度上解方程,大大地减小了计算时间。
多尺度有限元方法是由Babuska等提出的。该法在宏观尺度上进行网格剖分,然后通过在每个单元里求解细观尺度的方程(构造线性或者振荡的边界条件)来获得基函数。从而把细观尺度的信息反应到有限元法的基函数里,使宏观尺度的解包含了细观尺度的信息。但多尺度有限元方法在构造基函数时需要较大的计算量。
以上介绍了计算机数值计算常用的几种方法,了解这些方法原理及其优缺点有助于我们加深对数值计算的认识,方便我们在数值计算过程中掌握主动权。