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

《有限差分法求解偏微分方程》由会员分享，可在线阅读，更多相关《有限差分法求解偏微分方程（16页珍藏版）》请在装配图网上搜索。

1、有限差分法求解偏微分方程摘要：本文重要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基本，并在此基本上给出了部分有限差分法求解偏微分方程的算例验证了推导的对的性及操作可行性。核心词：计算力学，偏微分方程，有限差分法Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived

2、 in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method.Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Meth

3、od1 引言机械系统设计常常需要从力学观点进行构造设计以及构造分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有有关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。求解这些数学模型的措施大体分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载状况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，重要分为差分法和积分法两种，本论文重要讨论差分法。2 有限差分法理论基本2.1 有限差分法的基本

4、思想当系统的数学模型建立后，我们面对的重要问题就是微分积分方程的求解。基本思想是用离散的只具有限个未知量的差分方程组去近似地替代持续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似替代，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的环节重要有如下几步：n 区域离散，即把所给偏微分方程的求解区域细提成由有限个格点构成的网格，这些离散点称作网格的节点；n 近似替代，即采用有限差分公式替代每一种格点的导数；n 逼近求解，换而言之，这一过程可以看作是用一种插值多项式

5、及其微分来替代偏微分方程的解的过程。从原则上说，这种措施仍然可以达到任意满意的计算精度。由于方程的持续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应当收敛于精确解，但由于机器字节的限制，网格步长不也许也没有必要获得无限小，那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此，在运用有限差分法时，除了要保证精度外，还必须要保证其收敛性。2.2 系统微分方程的一般形式（1）由于大多数工程问题都是二维问题，因此得到的微分方程一般都是偏微分方程，对于一维问题得到的是常微分方程，解法与偏微分方程类似，故为了不是

6、一般性，这里只讨论偏微分方程。由于工程中高阶偏微分较少浮现，因此本文仅仅给出二阶偏微分方程的一般形式，对于高阶的偏微分，可进行类似地推广。二阶偏微分方程的一般形式如下：Axx+Bxy+Cyy=f(x,y,x,y)其中，为弹性体上的某一特性物理量（持续函数）。当A、B、C都是常数时，（1）式称为准线性，有三种准线性方程形式：n 如果=B2-4AC0，则称为双曲型方程。椭圆型方程重要用来解决稳态或静态问题，如热传导等问题；抛物线方程重要用来解决瞬态问题，如渗入、扩散等问题；双曲型方程重要用来解决振动问题，如玄震动、薄膜震动等问题。除了上述微分方程外，必须给出定解条件，一般有如下三类：n 第一类边界

7、条件（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=01n!fn(x0)x-x0n（2） =fx0+fxx-x0+12!fxx-x02+13!f(3)xx-x03+（3）对一种单变量函数fx

8、，以步长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-fxih+12!fxih2+R2(x)（4）fxi+h=fxi+fxih+12!fxih2+R2(x)由式（3）和式（4）可得到：（5）n 一阶向前差分：fxi=fxi+h-fxih（6）n 一阶向后差分：fxi=fxi-fxi-hh（7）n 一阶中心差分：fxi=fxi+h-fxi-h

9、2h不妨，记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-fih=fi+2-2fi+1+fih2（12）n 二阶向后差分：fi=fi-fi-1h=fi-2fi-1+fi-2h2（13）n 二阶中心差分：fi=fi+1-fi-12h=fi+2-2fi+fi-24h2差分公式（13）是以相隔2h的两结点处的函数值来表达中间结点处的一阶导数值，可称为中点导

10、数公式。式（11）和式（12）是以相邻三结点处的函数值来表达一种端点处的一阶导数值，可称为端点导数公式。应当指出：中点导数公式与端点导数公式相比，精度较高。由于前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽量应用前者，而只有在无法应用前者时才不得不应用后者。（14）但是，由于式（11）中的各阶导数均使用的是向前差分，导致用到的节点不相邻，同步为了均衡误差，将节点i处用到的一阶差分换成向后差分，则式（11）修正为：fi=fi+1-fih=fi+1-fih-fi-fi-1hh=fi+1-2fi+fi-1h2同理，根据上述推导过程，可得到任意阶的差分公式：（15）

11、n n阶向前差分：fi(n)=fi+1(n-1)-fi(n-1)h（16）n n阶向后差分：fi(n)=fi(n-1)-fi-1(n-1)h（17）n n阶中心差分：fi(n)=fi+1(n-1)-fi-1(n-1)2h阐明，上述公式中各节点处前一阶导数的代入也许存在不一致，也许是向前差分、向后差分或者中心差分，从而使最后的公式在系数上存在差别。固然，也可以对各相邻节点进行需要阶数的泰勒展开，从而建立方程组直接求各阶导数。2.3.2 微分方程转化为线性方程ym（18）由于三种类型的微分方程解法类似，故这里仅以椭圆型微分方程为例，将微分方程转化为代数方程，对于双曲型和抛物型方程依次类推即可。不妨

12、记：2u=uxx+uyy（2称为拉普拉斯算子），fx,y和g(x,y)是求解域上的持续函数。假设求解区域为：R=x,y:0xa,0yb,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-1y1 x1 xi-1 xi xi+1 xn图1 五点差分公式式（18）也称为五点差分公式，同理根据式

13、（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-2y

14、j-1yjyjy1yj+1yj+2 ui-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,jui+1,j+ui-1,j+ui,j+1+ui,j-1-4ui,jh2根据椭圆方程的具体形式可以将其分为如下三种形式：n 拉普拉斯（Lapl

15、ace）方程：2u=0n 泊松（Poison）方程：2u=g(x,y)n 赫耳墨次（Helmholtz）方程：2u+fx,yu=g(x,y)根据式（21），可建立三种不同形式椭圆方程的代数方程如下：n 拉普拉斯方程：2u=00=2ui,jui+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-h2gi,j=0n 赫耳墨次方程：2u+fx,yu=g(x,y)（24）ui

16、+1,j+ui-1,j+ui,j+1+ui,j-1+(h2fi,j-4)ui,j-h2gi,j=02.3.3 建立有限差分方程组根据式（22）（24）建立方程组，但是需要懂得相应的边界条件才干使方程组存在定解，根据2.2中可知，边界条件一般分为狄利克雷边界条件和导数边界条件两种，下面分别给出这两种边界条件的有限差分方程组的建立过程：n 狄利克雷边界条件：|=(x,y)对于狄氏条件而言，给出了边界上各节点出的函数计算公式，直接代入节点值(xi,yi)计算即可，如下所示为矩形区域的边界点计算：u(x1,yj)=u1,j=(x1,yj)(1jm)（左边界）（25）u(xn,yj)=un,j=(xn,

17、yj)(1jm)（右边界）u(xj,y1)=uj,1=(xj,y1)(1jn)（下边界）u(xj,yn)=uj,n=(xj,yn)(1jn)（上边界）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,yjun+1,j-un-1,j2

18、h根据式（28）可得到逼近表达：un+1,jun-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（左边界）图4 Neumann条件算子对于泊松方程和赫耳墨次方程同样根据上述措施，获得边界条件的线性方程，然后将这些方程添加到式（22）（24）所建立的方程组中，从而建立起(n-1)个(m-1)元的

19、线性方程组，解该方程组即可获得各节点的函数值。对于上述过程建立的线性方程组的求解，可采用多种措施，例如Jacob迭代法、Gauss-Seidel迭代法、超松弛迭代法（SOR法）、高斯消元法等措施求解。2.4 有限差分法的收敛性和稳定性由于迭代法必须保证收敛性，因此在解有限差分方程组时还应保证其收敛性，也就是一般所说的算法稳定性。有限差分法的算法稳定性可以通过特性值措施、傅里叶变换（冯诺依曼条件）以及能量估计等措施来判断，下面给出常用的冯诺依曼条件：n 向前差分：r1，绝对收敛；n 向后差分：r0，绝对收敛；n 中心差分：对任何的r对不收敛；假设求解域内x方向网格划分的步长为h，y方向网格划分的

20、步长为k，将偏微分方程化为原则形式，具体来说原则形式如下：（31）n 双曲方程：2uy2=c22ux2(c0)对于式（31）所示的双曲方程，冯诺依曼条件为：r=ckh.（32）n 抛物方程：uy=a2ux2对于式（32）所示的抛物方程，冯诺依曼条件为：r=ckh2.（33）n 椭圆方程：c22ux2+2uy2=0对于式（33）所示的椭圆方程，冯诺依曼条件为：r=ckh.（34）为了使算法在任何状况下都能保持稳定性，去掉对网格划分的冯诺依曼条件，一般采用隐式方案，对五点差分公式中的节点所在的行做差分，然后把这些差分的加权作为中心点的差分值，则拉普拉斯算子可修正为：uxx=ui+1,j+1-2ui

21、,j+1+ui-1,j+1h2+1-2ui+1,j-2ui,j+ui-1,jh2+ui+1,j-1-2ui,j-1+ui-1,j-1h2（01）运用式（34）进行计算时，稳定性没有任何限制。取不同的值得到不同的差分公式，一般取=1/4.（35）为了提高计算精度，很明显的一种措施就是网格细分，但是由于随着网格步长的减小，未知量的数目将会呈指数增长，网格划分太细会导致计算量过于庞大而无法计算。一般，我们可以通过提高逼近的精度，采用更高精度的差分公式，例如对公式（21）进行修改，可得到九点差分公式：2ui,j16h2ui+1,j+1+ui-1,j+1+ui+1,j-1+ui-1,j-1+4ui+1,

22、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 有限差分法求解实例根据上述推导有限差分法理论，对于不同类型的偏微分方程建立有限差分方程组，采用mat lab编程给出某些计算实例如下：1. 椭圆型方程n 拉普拉斯方程：2u=0；求解域：R=x,y:0x1.5,0y1.5下面分别给出拉普拉斯方程在不同的边界条件下的解。a) 狄利克雷边界条件：下边界：ux,0=x4上

23、边界：ux,1.5=x4-13.5x2+5.0625左边界：u0,y=y4右边界：u1.5,y=y4-13.5y2+5.0625图6 狄利克雷边界条件下拉普拉斯方程的解b) Neumann边界条件：u(x,y)N=0下边界：ux,0=0上边界：ux,1.5=400左边界：u0,y=0右边界：u1.5,y=0图7 Neumann边界条件下拉普拉斯方程的解n 泊松方程：2u=x2-y2；求解域：R=x,y:0x2,0y2.狄利克雷边界条件：下边界：ux,0=x4上边界：ux,1.5=x4-13.5x2+5.0625左边界：u0,y=y4右边界：u1.5,y=y4-13.5y2+5.0625图8 狄

24、利克雷边界条件下泊松方程的解n 赫耳墨兹方程：2u+(x+y)*u=x2-y2；求解域：R=x,y:0x1.5,0y1.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 狄利克雷边界条件下赫尔墨兹方程的解2. 双曲方程：uttx,t=4uxxx,t；求解区域：R=x,y:0x1,0t1.边界条件：u0,t=0u1,t=0ux,0=fx=sinx+sin2xuttx,0

25、=gx=0图10 双曲方程的解3. 抛物方程：utx,t=uxxx,t；求解区域：R=x,y:0x1,0t0.1.初值条件：ux,0=fx=sinx+sin3x边界条件：u0,t=g1t=0u1,t=g2t=0图11 抛物方程的解参照文献1. John H. Mathews, Kurtis D. Fink. Numerical Methods Using MATLAB, Fourth Edition. BEIJING: Publishing House of Electronic Industry, .72. 孙志忠. 偏微分方程数值解法（第二版）.北京：科学出版社，3. 王书亭.计算固体力学课件.