【《发展方程的新型数值计算方法及其应用概述》7700字】

发展方程的新型数值计算方法及其应用概述目录TOC\o"1-3"\h\u16161发展方程的新型数值计算方法及其应用概述 179581.1有限元法 1164861.1.1Galerkin有限元法 2294051.1.2修正的Galerkin有限元法 5190601.1.3特征有限元法 782371.2全离散计算格式 9197701.2.1简单的全离散格式 9123751.2.2高阶精度单步格式 10105891.2.3质量集中方法 11296061.3有限体积法 13244261.4谱与拟谱法 16154101.4.1投影与差值算子的逼近性质 1614911.4.2谱方法与拟谱方法 1821651.5保结构算法 221.1有限元法二阶抛物方程用于描写热的传播、溶质在液体中的弥散、多孔介质中渗流等随时间发展(演化)的现象和过程，是一类基本的发展型偏微分方程.求解抛物方程的初边值问题在科学与工程中有着广泛的应用。有限差分法可用于并曾经是近似求解抛物问题的数值方法，但在空间区域几何形式或者边值条件比较复杂的情形，该方法显得很笨重并且难于提高方法的精确度。由于抛物方程的初边值问题与椭圆方程的边值问题之间的密切关系，在近代有限元法诞生并成功的用于求解椭圆方程的边值问题之后，不久便被推广应用于求解抛物方程的初边值问题，并建立起完善的数学理论。所以我们以二阶抛物方程为例，介绍求解二阶线性抛物问题的有限元方法。设Ω是d维欧氏空间Rd中的一个有界区域，其边界记为Γ。考虑如下二阶线性抛物方程(1.1)其中A是由下式给出的一个微分算子如果下述条件成立，即α为某一正常数，对一切x∈Ω和ξ∈Rd，则称A为一致强椭圆微分算子，此时(1.1)属于抛物型方程。其次给定如下Dirichlet型边值条件(齐次)(1.2)和初值条件(1.3)1.1.1Galerkin有限元法本节讨论抛物方程初边值问题的Galerkin有限元近似（空间变量的离散化）。首先，借助虚功原理将问题改成变分形式。记设f∈C(0,T;H)，用函数v(x)∈H01与方程(利用Green公式因v=0当x∈Γ，故上式右端第二项为零。引进双线性泛函那么可得如下方程：(1.4)我们称方程(1.4)为问题(1.1)～(1.3)的变分(弱)形式，显然，初边值问题(1.1)～(1.3)的解u必定满足(1.4)，而问题(1.4)的解u(t):[0,T]→H01则称为原问题(1.1)~(1.3)的弱解，它被视作一个从I=[0,T]到空间由假设:aij=aji和定义式(1.6),对任意u,v∈H1,a(u,v)=a(v,u)，即a(·,·)是对称的双线性泛函。另外，根据关于算子A是一致强椭圆微分算子的假设，可知存在常数γ>0,使得即a(·,·)在空间H01上是正定的.从a(·,·)的定义还可看出，当aij(x),c(x)∈C(其中M为某一正常数.而有限元近似的第一步，是将求解空间区域Ω剖分成有限个互不重叠的子区域(称为元素或单元)的集合，用h表示剖分中单元的最大直径，记相应的剖分为Jh，{Jh,0<h<h0}代表一个剖分族，当h减小时，剖分就愈来愈细.以ρ(K)和σ(K)分别表示Jh中单元K的外接圆和内切圆的直径，如果存在不依赖于h的常数C,使得则称剖分族{Jh,0<h<h0}是正则的。有限元近似的第二步，是构造时的一族有限维子空间{Sh,0<h<h0},要求它具有如下逼近性质:对于某一整数τ≥2,有(1.5)通常，Sh是通过在剖分Jh上作分片多项式插值的方法去构造的，由Sobolev空间插值理论[6]可知，当剖分族{Jh,0<h<h0}为正则时，由Jh上所有属于C(Ω)的、次数≤r-1的分片多项式组成的子空间Sh满足条件(1.5)。在取定有限元空间Sh⊂H01后，初边值问题(1.1)~(1.3)的有限元近似定义为:求映射uh(t):I=[0,T]→(1.6)其中uhO∈Sh是函数u0(x)的某个近似。这里，我们看到问题(1.6)是变分问题(1.设{ϕj(x)}j=1N即确定其中系数{α(1.7)其中αj=dαjdt,rj为uh0(x)=j=1Nhrjϕ引进矩阵记号：其中M是一个Gram矩阵，称为“质量矩阵”，它是非奇异的，故M-1存在。利用上述记号，半离散问题(1.6)又可以写成(1.8)由常微理论可知，处置问题(1.7)对于任意F和r存在惟一的解α(t)，从而近似问题(1.6)存在唯一解uh(x，t).最后补充一点，即当初值问题中的边值条件不是Dirichlet型条件时，那么需要适当的修改变分陈述。例如，边值条件属于第三类条件：的情形，变分陈述(1.4)中的a(·,·)应当修改为并且尚需将其中的空间H01改换为1.1.2修正的Galerkin有限元法前面讨论的Galerkin有限元法，要求有限元近似空间Sh⊂H01(Ω),即Sh中的函数必需满足边值条件：在Γ=∂Ω上恒等于零，当Ω为带曲边界的区域时，构造满足上述要求的空间Sh是比较难的。为了克服这一困难，人们提出了种种修正的Galerkin有限元法，如拉氏乘子法[7]、插值边界条件法[8]和Nitsche方法[9]等。这些方法的共同特点是仅仅要求近似空间Sh⊂本节考虑如下含低阶微商项的二阶线性抛物方程的初边值问题：其中拉氏乘子法此方法是Babuska建议的。设和分别是H1(Ω)和H1(∂Ω)的子空间和逆性质用<·,·>代表L2(∂Ω)中的内积。令δ同时，定义椭圆近似算子Th：L2(Ω)对于上述{Th,0<h<h0}需要补充证明其中nj是∂Ω的单位外法向量的第j个分量。由迹定理其次，根据Sh的定义以及和的逼近性质，适当的选取，可得证毕。Nitsche方法为解除边界上限制，此方法利用了一个修正的双线性泛函代替a(v,w)，其中β为某一取定的正常数。引进范数假定子空间Sh具有逼近性质和如下迹不等式可以证明，对适当选定的常数β，ah(v,w)满足：借助Sh和双线性泛函ah(·,·)，椭圆近似算子Th：L2(Ω)→这类似于建立如下的误差估计：其次，由Th的定义，有其中等号成立仅当Thf=0。若vh∈Sh且(vh,Thvh)=0，则故vh=0。下令vh=Thf，wh=Thg，则有由|||·|||的定义和迹不等式，对于上式右端第二项，有很明显，对于右端的第一项也有这样的上方估计。1.1.3特征有限元法基于对流占优扩散问题具有双曲方程的特点，而存在特征方向(线)是双曲方程的重要特性。因此，在双曲方程和对流占优扩散方程的数值解法研究中，有许多工作致力于将普通的差分方法、有限元方法与特征线法结合起来见[10]。下面介绍特征线法和有限元方法的结合。以线性对流-扩散问题为例(1.9)在空间(x1,x2,t)中，由(1.10)所确定的曲线族称为方程(1.9)的特征(严格讲应属于μ=0的方程)。令ψ=(1+b12+借助特征方向，可将方程(1.9)改写为(1.11)首先，构造方程(1.11)的关于时间变量的差分逼近.令△t>0为时间步长，记tn=n△t,n=0,1,…由(1.10)知，t=tn时刻从x=(x1,x2)点始发的特征我在空间(t,x1,x2)与平面t=tn-1的交点为因此，这里用如下近似公式逼近沿特征方向的导数记u(x,tn)=un(x),由此可建立方程(1.9)的下述(时间离散化)近似方程(1.12)下一步要作的事情是，采用Galerkin有限元法，将(1.12)关于空间变量x进一步离散化。设{Jh,0<h<h0}为区域S的正则剖分族，Sh⊂H01(Ω)是基于剖分Jh和分片多项式插值构造的有限元函数空间，并假定Sh利用方程(1.12)的变分(弱)形式和子空间Sh。求解对流-扩散方程(1.9)的特征有限元法定义为:求函数uhn(x)∈Sh，n=0,1…1.2全离散计算格式抛物型偏微分方程的初边值问题，经过前文介绍的有限元法近似转化为常微分方程组的初值问题，其中时间变数仍为连续变量，所以称为半离散问题。为了获得所期待的数值解，尚需将前章所得半离散问题对时间变量进一步离散化，即建立全离散计算格式。半离散问题的完全离散化可以通过差分逼近、指数函数的有理逼近等方法来完成。1.2.1简单的全离散格式考虑热传导问题：(1.13)其中为laplace算子。在取定有限元近似空间Sh⊂H01(Ω)之后，上式的Galerkin有限元(半离散)近似为：求映射uh(t):I=[0,T]→(1.14)本节采用差分逼近方法构造问题(1.13)的两个常用的全离散格式。为此，用节点0=t0<t1<…<tN=T代替时间区间I=[0,T],记In=(tn-l,tn)，称kn=tn-tn-1为时间步长，目的是求初边值问题真解u(x,t)在节点tn,n=1,2,…处的近似值，记作Un.向后Euler-Galerkin格式求Un∈Sh，n=0,1,2,…使得(1.15)此格式是由向后差商(un(tn)-uh(tn-1)/kn近似方程(1.14)中的导数∂uh(tn)/∂t,再将uh(tl)改记为Ul而导出的。不难看出，(1.15)相对于方程(1.14)的局部截断误差为O(kn).采用1.1中的矩阵记号，格式(1.15)可改写成(1.16)由此看到，当Un-1(即αn-1)已知，那么通过求解具正定、对称系数矩阵(M+knB)的线性代数方程组(1.16)，便可求出下一时刻的近似解Un,这是一个隐型格式。Crank-Nicolson-Galerkin格式求Un∈Sh，n=0,1,2,…，使得此格式是由在tn−12=tn−12kn点用中心差商(uh(tn)-uh(tn-1))/kn近似(1.14)中的导数∂uh(tn−12)/∂t并用平均值(uh(tn)-uh(tn-1))/2近似采用矩阵记号，此格式可以改写为这也是一个隐格式。1.2.2高阶精度单步格式本节介绍通过指数函数的有理逼近来建立半离散问题的高阶精度单步格式的方法。先考虑齐次热传导问题：设{Sh,h>0}为有限维空间族，并给定了逼近T=(-∆)-1的一族算子{Th,h>0}，其中算子Th:L2(Ω)→STh是自共轭的，它在L2(Ω)上半正定并在S对任意的f∈Hs(Ω)，0≤Th代表一般的椭圆逼近算子。然后令△h=-Th-1，并视△h为从Sh到Sh的算子(离散Laplace算子)。由此，可以定义上式的半离散近似：求uh(t)∈Sh,t≥0满足(1.17)此问题的解可表示为uh(t)=et△huh0现设r(z)是指数函数ez的一个有理分式近似：其中v为某一非负整数，并且r(z)在负实轴上无极点。利用函数r(z)，我们可以定义问题(1.17)的单步全离散格式：(1.18)其中Un,uh0均属于Sh，k>0为时间步长，并称(1.18)为(1.1.2.3质量集中方法本节介绍一个可以大大简化求解计算的修正Galerkin有限元法：质量集中方法。由于这个方法不仅使求解计算量显著降低，同时具有很好的稳定性和收敛性，所以受到人们的重视。回忆前文中中所定义的Galerkin有限元法，以热传导问题为例，在取定有限元空间Sh后，此问题的Galerkin有限元近似为(1.19)设{ϕj(x)}j=1N则问题(1.19)的矩阵形式为(1.20)其中M=(mij)Nh×Nh和B=(bij)Nh×Nh分别称为“质量矩阵”和“刚度矩阵”，它们的元素为mij=(ϕi;ϕj)和bij=(△ϕi,△ϕj)。若采用向后差商近似(2.43)中的微商dα(tn)/dt，可构造出如下全离散格式：(1.21)此即本章第一节中所述向后Euler-Galerkin格式，我们看到，在这里无论是(1.20)还是(1.21)的求解均要求计算逆矩阵M-1。在有限元法中，M通常是一个高阶的非奇异矩阵，计算M-1B和M-1F所需运算量和存储量是很大的。所谓的质量集中法，就是将(1.20)和(1.21)中的矩阵M(质量矩阵)用一个对角矩阵M=diag(m11,换言之，M是由将M中每行元素集中起来(作代数和)置于该行的对角线位置上。此时，方程(1.20)和(1.21)分别地简化成和显然，这样做可以大大简化求解计算。因为计算M-1要比计算M-1简单得多。设E(t,s)为热传导方程的解算子，即u(t)=E(t,s)us，t∈[s,+∞)代表齐次热传导方程满足初始条件：u|t=s=us的解。由抛物方程的极值性质，我们知道，算子E(t,s)满足由此可见，齐次热传导方程的解u(x,t)在最大值模的意义下是稳定的，即有质量集中方法的一个重要优点，就是当剖分满足一定条件时，相应的离散方程仍然具有上述极值性质。然而，对于通常的Galerkin有限元法，连续方程原有的这一性质会遭到破坏。1.3有限体积法我们以扩散问题为例介绍有限体积法。设Ω⊂R2是一个有界区域，x=(x1,x2)∈Ω。考虑如下线性对流-扩散问题：(1.22)其中μ>0和b1,b2为常数，并设μ<<|b|。设Ωh={K}为空间区域的正则三角剖分，K表示剖分Ωh中的三角形单元，Ωh中所有三角形顶点组成的节点集合为{Pi:1≤i≤M2}，其中{Pi:1≤i≤M;}为内节点集合，而{Pi:M1+1≤i≤M2}为边界节点集合。h代表Ωh中所有单元的最大直径。设函数空间Uh⊂H01(Ω)为剖分Ωh上的标准线性有限元空间，其基函数为{ϕi(x)}j=1M1满足条件ϕi(Pj)=δij,1≤i,j≤其次，按下述方法建立Sh的对偶剖分。对于任一三角形单元K∈Ωh,在K内选定一点Q，并作此点与K的各边中点的连线，则K被分成三个小三角形.对任意节点Pi,考虑所有以Pi为公共顶点的单元，将所有以Pi为顶点的小三角形的并集记为KPi∗,它是一个多边形子区域，被称为对偶单元。所有对偶单元的集合便构成Ωh的一个对偶剖分，记为Ωh∗，当Q点选作K的重心时，相应的Ωh∗称为重心对偶制分。当令Vh代表Ωh∗上所有分片常数且于KPi∗(M1+1≤i引进插值投影算子Πh:C0(Ω)→以及Πh∗:C0(Ω)→Ωh中两个节点Pi和Pj如果它属于某个单元边的端点，则称Pi和Pj为邻点.设Pi和Pj为相邻节点，记Γij=∂KPi∗∩∂K有限体积法的建立过程是这样的：首先，在KPi∗(称控制体积)上积分方程(利用Green公式v：∂KPi∗上的单位外法向量并将其中的u改为uh,得到如下积分形式的方程(1.23)由于un是Ωh上的分片线性函数，采用分片和分段积分的方式将上式左右两端的积分计算出来，给出它们的表达公式。如此，便可建立近似求解问题(1.22)的有限体积法的计算格式(半离散方程)。下面，就Ωh∗为外心对偶剖分的情形讨论一下(1.6)的具体形式，令Λi={j：Pj与Pi(1.24)令并将uh改记为uh(P,t)，利用上述积分表达式与近似公式，便可得到如下近似求解的计算格式另外，若引进Heaviside函数h(r)：并将前面关于对流项的近似公式(1.24)改为其余各项的近似保持不变，此时所得到的格式称为迎风型有限体积格式。出于理论分析和推广的需要，下面介绍方程(1.23)的变分形式。引进记号则可将(1.6)写成因此，求解方程(1.23)等价于：求uh(·,t)∈Uh满足:(1.25)此即方程(1.23)的变分形式，在变分问题(1.13)中，Uh称为试探函数空间(近似解所属空间)，而Vh则称为检验函数空间(用来检验解是否满足方程)。此类离散化称为Petrov-Galerkin有限元法，它与通常Galerkin有限元法的不同之处，就是其中的试探函数空间和检验函数空间是不相同的空间。1.4谱与拟谱法谱与拟谱法是以三角多项式(或其他正交多项式)为基函数的Galerkin方法或配置法，它在快速Fourier变换出现之后获得迅速的发展，成为近似求解偏微分方程的重要方法之一。1.4.1投影与差值算子的逼近性质Fourier(三角多项式)逼近设SN=span{eikx,-N≤k≤N-1},则SN为次数≤N的三角多项式所构成的函数空间，其维数是2N。用Lp2(0,2π)表示由2π对任意u∈Lp2(0,2π)，存在PNu∈称PNu为u的Fourier投影，算子PN:Lp2(0,2π)→SN称为Fourier投影算子。由函数系{eikx}的正交性:(eikx,eik’x则其中ϕk(x)=eikx，这表明Fourier投影PNu实际上就是uChebyshev多项式逼近作为Fouier逼近的一个重要补充，介绍Chebyshev多项式逼近，它对于求解非周期边值问题有特殊的用途[11]。Chebyshev多项式为此处，T0(x)=1,T1(x)=x,T2(x)=2x2-1,…。利用变换x=cosθ和公式可得如下递推关系式已知，{Tk(x)}是区间[-1,1]上以ω(x)=(1-x2)−12为权的正交函数系。事实上，若令x=因为故知：当k为偶数时，Tk(x)是偶函数，而当k为奇数时，Tk(x)是奇函数。另外，Tk(x)有下列性质：并有表达式：1.4.2谱方法与拟谱方法本节以热传导方程的初边值问题为例介绍谱与与拟谱方法，包括用Chebyshev配置法求解非周期边值条件的问题。这些方法都是用于空间变量的离散化，即构造问题的半离散近似。Pourier-Galerkin方法(谱方法)这里，考虑热传导方程的周期初值问题:(1.26)其中f(·,t),u0(·)∈Lp2(0,2π)。当u0∈Hps(0,2π)和f设SN=Span{eikx:-N≤k≤N-1},则问题(1.26)的Fourier-Calerkin近似为求uN(t)∈SN,使得(1.27)若令将它代人(1.27)并取v=eikx，由{eilx}l=−N它是一个以{Ck(t)}l=−NN−1为未知函数的一阶线性常微分方程组，求出系数{Ck(t)}l=−NN−1后，便可得到问题(1.Fouier配置法(拟谱方法)在区间[0,2π]上选取节点(配置点)xj=jh，j=0,1,…,2N-1,其中h=2π2N。则求解(1.26)的Fouier配置法为求uN(·,t)∈SN(1.28)此即要求近似解uN(x,t)在各个配置点处满足方程和初值条件。利用离散内积和范数：定义式(1.28)又可以改写为(1.29)若令其中(1.30)则有(1.31)和将uN和(uN)xx的上述表示式代入(1.29)，可得关系式(1.30)和(1.31)给出了{uN(xj,t)}j=02N−1(物理向量)和{uN,k(t)}K=−NN−1(谱向量)之间的转换公式，t被视作参数。通常称(1.30)为离散Fourier变换，而(5.39)则是它的逆变换，离散变换(比较一下定义式(1.28)和(1.29)，可以看出，拟谱方法与谱方法的不同之处就是将(1.28)中的内积(·,·)改为离散内积(·,·)h。Chebyshev配置法前面用Fourier-Galerkin法和Fouier配置法构造出含周期边值条件：u(x+2π,t)=u(x,t)的热传导(抛物方程)问题的离散近似。如何将谱(拟谱)方法推广，使之适用于非周期边值条件问题是人们颇为关注的研究课题。现在考虑如下形式的初边值问题:(1.32)其中Bu=0代表在区间[-1,1]端点x=±1处指定的边值条件(一般的属于非周期条件)。作为特例，Bu=u情形为Dirichlet边值条件，而情形Bu=ux则属Neumann边值条件。于此，在区间[-1,1]上选择配置点：xj=cosjπN，j=0,1,…,N.这是以ω然后，按配置法的原理定义初边值问题(1.32)的下述离散化近似：求uN(·,t)∈PN,t∈[0,T]，使得(1.33)设其中下面讨论uN(x,t)对x微商的Chebyshev展开。为简化记号，记则因u(1)(x)∈PN-1，故可令由可知因uk(1)=0一般的，有特别由可得二阶导数u(2)(x)的Chebyshev展开。将uN(xj,t)及其导数的展开式代入(1.33)，则得关于{uN,k(t)}K=0其中(u0)k中的系数，当解此初值问题求出{uN,k(t)}K=0即可得所需的近似解，于此可以借助类似于离散Fouier变幻的快速算法得以实现。最后，仅就Bu=u即Dirichlet边值条件：u(±1,t)=0的情形。1.5保结构算法动力系统的几何算法诞生于20世纪80年代中期[12]，冯康从哈密顿系统相空间的辛几何结构出发提出了辛算法(格式)，并在国外倡导了几何算法的研究。发展型偏微分方程可视为无穷维的动力系统，人们已发现许多物理、力学中的重要非线性发展方程同样地具有辛几何结构(或Lie群结构)，例如Schrodinger方程、Sine-Gordon方程、Boussinisq和KortewegdeVries方程等。在几何算法的思想启示下，近十余年来出现了一批关于非线性发展方程的保几何结构算法的研究工作，目前已成为发展方程数值方法研究的热点之一。下以几类重要哈密顿型非线性发展方程介绍如何构造它们的保结构数值方法。在经典力学中，含n个自由度的，受位势力作用的系统(称保守系统)均可用如下形式的一阶常微分方程组描述：(1.34)其中H=H(q1,…,qn,p1,…,pn)称为哈密领函数。它代表所描述力学系统的总能量，q(t)=(q1(t),…,qn(t))为系统的广义坐标。p(t)=(p1(t),…,pn(t))为广义动量。方程组(1.34)称为n自由度力学系统的哈密顿形式，数学上称它为n自由度哈密顿系统。描述n自由度系统的运动的熟知形式为牛顿形式和拉格朗日变分形式。对于保守系统，设势能函数为V=V(q)，则方程为保守系统的牛顿形式。另外，设L=L(q,q,t)(=动能