有限差分法求解偏微分方程

《有限差分法求解偏微分方程》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

有限差分法求解偏微分方程摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。关键词：计算力学，偏微分方程，有限差分法Abstract：Thisdissertationmainlyfocusesonsolvingthemathematicmodelofcomputationmechanicswithfinite-differencemethod.Thetheoreticalbasisoffinite-differenceisderivedinthesecondpartofthedissertation,andthenIuseMATLABtoprogramthealgorithmstosolvesomepartialdifferentialequationstoconfirmthecorrectnessofthederivationandthefeasibilityofthemethod.Keywords：ComputationMechanics,PartialDifferentialEquations,Finite-DifferenceMethod 1引言机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，主要分为差分法和积分法两种，本论文主要讨论差分法。2有限差分法理论基础2.1有限差分法的基本思想当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的步骤主要有以下几步：n区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格，这些离散点称作网格的节点；n近似替代，即采用有限差分公式替代每一个格点的导数；n逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小， 那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此，在运用有限差分法时，除了要保证精度外，还必须要保证其收敛性。2.2系统微分方程的一般形式（1）由于大多数工程问题都是二维问题，所以得到的微分方程一般都是偏微分方程，对于一维问题得到的是常微分方程，解法与偏微分方程类似，故为了不是一般性，这里只讨论偏微分方程。由于工程中高阶偏微分较少出现，所以本文仅仅给出二阶偏微分方程的一般形式，对于高阶的偏微分，可进行类似地推广。二阶偏微分方程的一般形式如下：Aϕxx+Bϕxy+Cϕyy=f(x,y,ϕ,ϕx,ϕy)其中，ϕ为弹性体上的某一特征物理量（连续函数）。当A、B、C都是常数时，（1）式称为准线性，有三种准线性方程形式：n如果Δ=B2-4AC<0，则称为椭圆型方程；n如果Δ=B2-4AC=0，则称为抛物型方程；n如果Δ=B2-4AC>0，则称为双曲型方程。椭圆型方程主要用来处理稳态或静态问题，如热传导等问题；抛物线方程主要用来处理瞬态问题，如渗透、扩散等问题；双曲型方程主要用来处理振动问题，如玄震动、薄膜震动等问题。除了上述微分方程外，必须给出定解条件，通常有如下三类：n第一类边界条件（Dirichlet条件）：ϕ|Γ=φ(x,y)；n第二类边界条件（Neumann条件）：∂ϕ∂n|Γ=φ1(x,y)；n第三类边界条件（Robin条件）：[∂ϕ∂n+λ(x,y)ϕ]|Γ=ψ(x,y)；其中，Γ为求解域Ω的边界，n为Γ的单位外法矢，λ(x,y)|Γ≢0。第二类和第三类边界条件统称为导数边界条件。2.3有限差分方程的数学基础2.3.1一元函数导数的差分公式一个函数在x点上的导数，可以近似地用它所临近的两点上的函数值的差分 来表示。函数fx在x=x0处的泰勒展式如下：fx=n=0∞1n!fn(x0)x-x0n（2）=fx0+f'xx-x0+12!f''xx-x02+13!f(3)xx-x03+⋯（3）对一个单变量函数fx，以步长∆x=h将[a,b]区间等距划分，我们得到一系列节点：x0=a，x1=x0+∆x=a+h，x2=x1+∆x=a+2h，⋯，xi=xi-1+∆x=a+ih，⋯，xn=b（n=b-ah），然后求出fx在这些节点上的近似值。与节点xi相邻的节点有xi-h和xi+h，因此在点xi处可以构造如下形式的展开式：fxi-h=fxi-f'xih+12!f''xih2+R2(x)（4）fxi+h=fxi+f'xih+12!f''xih2+R2(x)由式（3）和式（4）可得到：（5）n一阶向前差分：f'xi=fxi+h-fxih（6）n一阶向后差分：f'xi=fxi-fxi-hh（7）n一阶中心差分：f'xi=fxi+h-fxi-h2h不妨，记fi=f(xi)，则式（5）、（6）、（7）分别简写为：（8）n一阶向前差分：fi'=fi+1-fih（9）n一阶向后差分：fi'=fi-fi-1h（10）n一阶中心差分：fi'=fi+1-fi-12h根据式（8）、式（9）和式（10），可得二阶差分： （11）n二阶向前差分：fi''=fi+1'-fi'h=fi+2-2fi+1+fih2（12）n二阶向后差分：fi''=fi'-fi-1'h=fi-2fi-1+fi-2h2（13）n二阶中心差分：fi''=fi+1'-fi-1'2h=fi+2-2fi+fi-24h2差分公式（13）是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。式（11）和式（12）是以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。（14）但是，由于式（11）中的各阶导数均使用的是向前差分，导致用到的节点不相邻，同时为了均衡误差，将节点i处用到的一阶差分换成向后差分，则式（11）修正为：fi''=fi+1'-fi'h=fi+1-fih-fi-fi-1hh=fi+1-2fi+fi-1h2同理，根据上述推导过程，可得到任意阶的差分公式：（15）nn阶向前差分：fi(n)=fi+1(n-1)-fi(n-1)h（16）nn阶向后差分：fi(n)=fi(n-1)-fi-1(n-1)h（17）nn阶中心差分：fi(n)=fi+1(n-1)-fi-1(n-1)2h说明，上述公式中各节点处前一阶导数的代入可能存在不一致，可能是向前差分、向后差分或者中心差分，从而使最终的公式在系数上存在差别。当然，也可以对各相邻节点进行需要阶数的泰勒展开，从而建立方程组直接求各阶导数。 2.3.2微分方程转化为线性方程ym（18）由于三种类型的微分方程解法类似，故这里仅以椭圆型微分方程为例，将微分方程转化为代数方程，对于双曲型和抛物型方程依次类推即可。不妨记：∇2u=uxx+uyy（∇2称为拉普拉斯算子），fx,y和g(x,y)是求解域上的连续函数。假设求解区域为：R={x,y:0≤x≤a,0≤y≤b,ba=m/n}，将求解区域划分成(n-1)×(m-1)个网格，其中：a=nh,b=mh，如图1所示。记fi,j=f(xi,yj)，则根据式（14）可得到：∇2u=uxx+uyy=ui+1,j-2ui,j+ui-1,jh2+ui,j+1-2ui,j+ui,j-1h2=ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2+O(h2)yj+1yjyj-1y1x1xi-1xixi+1xn图1五点差分公式式（18）也称为五点差分公式，同理根据式（12）和式（13）可分别得到向前差分公式（19）和向后差分公式（20），如图（2所示）。（19）n向前差分∇2u=uxx+uyy=ui+2,j-2ui+1,j+ui,jh2+ui,j+2-2ui,j+1+ui,jh2=ui+2,j+ui,j+2+2ui,j-2ui+1,j-2ui,j+1h2+O(h2)ymym（20）n向后差分 ∇2u=uxx+uyy=ui,j-2ui-1,j+ui-2,jh2+ui,j-2ui,j-1+ui,j-2h2=ui-2,j+ui,j-2+2ui,j-2ui-1,j-2ui,j-1h2+Oh2x1xi-2xi-1xixi+1xi+2xixnxnx1y1yj-2yj-1yjyjy1yj+1yj+2ui-2,j2ui,jui,j+2ui,j+1图2向前差分（左）和向后差分（右）-2ui-1,j-2ui,j-1ui,j-2-2ui,j+12ui,j-2ui+1,jui+2,j-4ui,jui,j-1ui+1,jui-1,j图3中心差分、向前差分和向后差分的拉普拉斯算子表示（21）利用中心差分公式（18），由于式（18）在点x,y=(xi,yi)处具有二阶精度（Oh2），所以式（18）可近似改写成下式：∇2ui,j≈ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2根据椭圆方程的具体形式可以将其分为以下三种形式：n拉普拉斯（Laplace）方程：∇2u=0n泊松（Poison）方程：∇2u=g(x,y)n赫耳墨次（Helmholtz）方程：∇2u+fx,yu=g(x,y)根据式（21），可建立三种不同形式椭圆方程的代数方程如下：n拉普拉斯方程：∇2u=00=∇2ui,j≈ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2化简后得到拉普拉斯方程的计算公式：（22）ui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,j=0 （23）n泊松方程：∇2u=gx,yui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,j-h2∙gi,j=0n赫耳墨次方程：∇2u+fx,yu=g(x,y)（24）ui+1,j+ui-1,j+ui,j+1+ui,j-1+(h2∙fi,j-4)ui,j-h2∙gi,j=02.3.3建立有限差分方程组根据式（22）~（24）建立方程组，但是需要知道对应的边界条件才能使方程组存在定解，根据2.2中可知，边界条件一般分为狄利克雷边界条件和导数边界条件两种，下面分别给出这两种边界条件的有限差分方程组的建立过程：n狄利克雷边界条件：ϕ|Γ=φ(x,y)对于狄氏条件而言，给出了边界上各节点出的函数计算公式，直接代入节点值(xi,yi)计算即可，如下所示为矩形区域的边界点计算：u(x1,yj)=u1,j=φ(x1,yj)(1≤j≤m)（左边界）（25）u(xn,yj)=un,j=φ(xn,yj)(1≤j≤m)（右边界）u(xj,y1)=uj,1=φ(xj,y1)(1≤j≤n)（下边界）u(xj,yn)=uj,n=φ(xj,yn)(1≤j≤n)（上边界）n导数边界条件：∂u(x,y)∂N=0（26）以右边界点为例，对于右边界点x=xn，根据Neumann条件可得下式：∂u(x,y)∂N=∂u(xn,yj)∂x=uxxn,yj=0对于拉普拉斯方程，根据计算公式（22），对于边界上的点(xn,yj)可得：（27）un+1,j+un-1,j+un,j+1+un,j-1-4un,j=0（28）显然，上式中的un+1,j在求解域外，是未知量。根据中心差分公式（10）可得到：uxxn,yj≈un+1,j-un-1,j2h根据式（28）可得到逼近表示：un+1,j≈un-1,j，并且具有2阶逼近精度，代入式（27）可得下式：（29）2un-1,j+un,j+1+un,j-1-4un,j=0 同理，对于其它边界可获得如下边界方程：（30）2ui,2+ui-1,1+ui+1,1-4ui,1=0（下边界）2ui,m-1+ui-1,m+ui+1,m-4ui,m=0（上边界）2u2,j+u1,j+1+u1,j-1-4u1,j=0（左边界）图4Neumann条件算子对于泊松方程和赫耳墨次方程同样根据上述方法，获得边界条件的线性方程，然后将这些方程添加到式（22）~（24）所建立的方程组中，从而建立起(n-1)个(m-1)元的线性方程组，解该方程组即可获得各节点的函数值。对于上述过程建立的线性方程组的求解，可采用多种方法，比如Jacob迭代法、Gauss-Seidel迭代法、超松弛迭代法（SOR法）、高斯消元法等方法求解。2.4有限差分法的收敛性和稳定性由于迭代法必须保证收敛性，所以在解有限差分方程组时还应保证其收敛性，也就是通常所说的算法稳定性。有限差分法的算法稳定性可以通过特征值方法、傅里叶变换（冯诺依曼条件）以及能量估计等方法来判断，下面给出常用的冯诺依曼条件：n向前差分：r≤1，绝对收敛；n向后差分：r>0，绝对收敛； n中心差分：对任何的r对不收敛；假设求解域内x方向网格划分的步长为h，y方向网格划分的步长为k，将偏微分方程化为标准形式，具体来说标准形式如下：（31）n双曲方程：∂2u∂y2=c2∂2u∂x2(c>0)对于式（31）所示的双曲方程，冯诺依曼条件为：r=ckh.（32）n抛物方程：∂u∂y=a∂2u∂x2对于式（32）所示的抛物方程，冯诺依曼条件为：r=ckh2.（33）n椭圆方程：c2∂2u∂x2+∂2u∂y2=0对于式（33）所示的椭圆方程，冯诺依曼条件为：r=ckh.（34）为了使算法在任何情况下都能保持稳定性，去掉对网格划分的冯诺依曼条件，通常采用隐式方案，对五点差分公式中的节点所在的行做差分，然后把这些差分的加权作为中心点的差分值，则拉普拉斯算子可修正为：uxx=θui+1,j+1-2ui,j+1+ui-1,j+1h2+1-2θui+1,j-2ui,j+ui-1,jh2+θui+1,j-1-2ui,j-1+ui-1,j-1h2（0≤θ≤1）利用式（34）进行计算时，稳定性没有任何限制。θ取不同的值得到不同的差分公式，通常取θ=1/4.（35）为了提高计算精度，很明显的一个措施就是网格细分，但是由于随着网格步长的减小，未知量的数目将会呈指数增长，网格划分太细会导致计算量过于庞大而无法计算。通常，我们可以通过提高逼近的精度，采用更高精度的差分公式，例如对公式（21）进行修改，可得到九点差分公式：∇2ui,j≈16h2ui+1,j+1+ui-1,j+1+ui+1,j-1+ui-1,j-1+4ui+1,j+4ui-1,j+4ui,j-1+4ui,j+1-20ui,j+O(h4) ui-1,j+14ui-1,jui-1,j-14ui+1,j-20ui,j4ui-1,j4ui,j+1ui-1,j+1ui+1,j+1xi+1xixi-1x1yj+1yj-1yjy1ym图5九点差分公式3有限差分法求解实例根据上述推导有限差分法理论，对于不同类型的偏微分方程建立有限差分方程组，采用matlab编程给出一些计算实例如下：1.椭圆型方程n拉普拉斯方程：∇2u=0；求解域：R={x,y:0≤x≤1.5,0≤y≤1.5}下面分别给出拉普拉斯方程在不同的边界条件下的解。a)狄利克雷边界条件：下边界：ux,0=x4上边界：ux,1.5=x4-13.5x2+5.0625左边界：u0,y=y4右边界：u1.5,y=y4-13.5y2+5.0625 图6狄利克雷边界条件下拉普拉斯方程的解a)Neumann边界条件：∂u(x,y)∂N=0下边界：ux,0=0上边界：ux,1.5=400左边界：u0,y=0右边界：u1.5,y=0图7Neumann边界条件下拉普拉斯方程的解 n泊松方程：∇2u=x2-y2；求解域：R={x,y:0≤x≤2,0≤y≤2}.狄利克雷边界条件：下边界：ux,0=x4上边界：ux,1.5=x4-13.5x2+5.0625左边界：u0,y=y4右边界：u1.5,y=y4-13.5y2+5.0625图8狄利克雷边界条件下泊松方程的解n赫耳墨兹方程：∇2u+(x+y)*u=x2-y2；求解域：R={x,y:0≤x≤1.5,0≤y≤1.5}.狄利克雷边界条件：下边界：ux,0=x4上边界：ux,1.5=x4-13.5x2+5.0625左边界：u0,y=y4右边界：u1.5,y=y4-13.5y2+5.0625通过图9与图6的对比发现，微分方程的解与微分方程的具体形式关系不大，主要由求解域和边界条件所决定。 图9狄利克雷边界条件下赫尔墨兹方程的解1.双曲方程：uttx,t=4uxxx,t；求解区域：R={x,y:0≤x≤1,0≤t≤1}.边界条件：u0,t=0u1,t=0ux,0=fx=sinπx+sin2πxuttx,0=gx=0图10双曲方程的解 1.抛物方程：utx,t=uxxx,t；求解区域：R={x,y:0≤x≤1,0≤t≤0.1}.初值条件：ux,0=fx=sinπx+sin3πx边界条件：u0,t=g1t=0u1,t=g2t=0图11抛物方程的解参考文献1.JohnH.Mathews,KurtisD.Fink.NumericalMethodsUsingMATLAB,FourthEdition.BEIJING:PublishingHouseofElectronicIndustry,2005.72.孙志忠.偏微分方程数值解法（第二版）.北京：科学出版社，20123.王书亭.计算固体力学课件.