偏微分方程数值解法

1、十 章 偏 微 分方 程 数 值 解 法偏微分方程问题，其求解十分困难。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 1差分方法的基本概念1.1几类偏微分方程的定解问题椭圆型方程：其最典型、最简单的形式是泊松( Poisson )方程特别地，当 f(x,y) 0时，即为拉普拉斯(Laplace )方程，又称为调和方程Poisson方程的第一边值问题为其中 为以 为边界的有界区域，为分段光滑曲线，称为定解区域，f (x, y), (x, y)分别为 ,上的已知连续函数。第二类和第三类边界条件可统一表示为其中n为

2、边界的外法线方向。当0时为第二类边界条件，0时为第三类边界条件。抛物型方程：其最简单的形式为一维热传导方程方程可以有两种不同类型的定解问题：初值问题初边值问题其中 (x), g1(t), g2(t)为已知函数，且满足连接条件边界条件 u(0,t) g1(t),u(l,t) g2(t)称为第一类边界条 件。第二类和第三类边界条件为其中 1(t) 0, 2(t) 0。当 1(t)2(t)0 时，为第二类边界条件，否则称为第三类边界条件。双曲型方程：最简单形式为一阶双曲型方程物理中常见的一维振动与波动问题可用二阶波动方程描述，它是双曲型方程的典型形式。方程的初值问题为边界条件一般也有三类，最简单的初

3、边值问题为1.2差分方法的基本概念差分方法又称为有限差分方法或网格法，是求偏微分方程定 解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连 续变化区域用有限离散点(网格点)集代替；将问题中出现的连 续变量的函数用定义在网格点上离散变量的函数代替；通过用网 格点上函数的差商代替导数，将含连续变量的偏微分方程定解问 题化成只含有限个未知数的代数方程组(称为差分格式)。如果 差分格式有解，且当网格无限变小时其解收敛于原微分方程定解 问题的解，则差分格式的解就作为原问题的近似解(数值解)。 因此，用差分方法求偏微分方程定解问题一般需要解决以下问题:(1)选取网格

4、；(2)对微分方程及定解条件选择差分近似，列出差分格式；(3)求解差分格式；(4)讨论差分格式解对于微分方程解的收敛性及误差估计。下面，用一个简单的例子来说明用差分方法求解偏微分方程 问题的一般过程及差分方法的基本概念。设有一阶双曲型方程初值问题。(1)选取网格：-2h-h0h2h3h首先对定解区域D (x,t) x ,t 0作网格剖分，最简单常用一种网格是用两族分别平行于 x轴与t 轴的等距直线x Xk kh,t tj j (k 0, 1, 2L ,j 0,1,2,L u等D 分成许 多小矩形区域。这些直线称为网格线，其交点称为网格点，也称为节点,h和 分别称作 X 方向和t方向的步长。这种

5、网格称为矩形网格(2)对微分方程及定解条件选择差分近似，列出差分格式： 如果用向前差商表示一阶偏导数，即其中01, 21方程处可表小为其中 u(Xk,0)(Xk)(k 0,1,2,L)。由于当h, 足够小时，在式中略去R(Xk,tj),就得到一个与方程相近似的差分方程此处，Uk, j可看作是问题的解在节点(x k , t j )处的近似值。同初值条件结合，就得到求问题的数值解的差分格式。式O( qh p),则称差分方称为差分方程的截断误差。如果一个差分方程的截断误差为 R程对t是q阶精度，对x是p阶精度的。显然，截断误差的阶数 越大，差分方程对微分方程的逼近越好。若网格步长趋于0时，差分方程的

6、截断误差也趋于0,则称 差分方程与相应的微分方程是相容u 这是用差分方法求解偏微 分方程问题的必要条件。如果当网格步长趋于0时，差分格式的解收敛到相应微分方 程定解问题的解，则称这种差分格式是收敛的。 2椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法。差分格式的建立考虑Poisson方程的第一边值问题取h, 分别为x方向和y方向的步长，如图所示，以两族平行线x xk kh, y yj j (k,j 0, 1, 2,L )将定R (xk, yj)解区域剖分成矩形网 格。节点的全体记为xk kh, yjj , k, j为整数 0定解区域内部的节点称为内

7、点，记内点集 R 为 h 。边界 与网格 线的交点称为边界点，边界点全体记为h o与节点(x, yj)沿x方向或y方向只差一个步长的点(xk i,yjM(xk,yj .称为节点(Xk,yj)的相邻节点。如果一个内点的四个相邻节点均属于，称为正一 (i)一 一则内点，正内点的全体记为，至少有一个相邻干点不属于-一“一 “一，(2)、目工的内点称为非正则内点，非正则内点的全体记为。问题是要求出第一边值问题在全体内点上的数值解。为简便，记(k,j) (Xk,yj), u(k, j) u(4,yj),fk,j f(xk，yj)。对正则(1),人，,内点(k,j),由二阶中心差冏公式2u2Poisson

8、 方程 x2了 f(x，V)在点(k, j)处可表小为 其中为其截断误差表示式，略去R(k, j),即得与方程相近似的差分方程式中方程的个数等于正则内点的个数，而未知数Uk,j则除了包含正则内点处解 U 的近似值外，还包含一些非正则内点处U的近似值，因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按上式给出，需要利用边界条件得到。边界条件的处理可以有各种方案，下面介绍较简单的两种。 (1)直接转移用最接近非正则内点的边界点上的 u值作为该点上 u 值的 近似，这就是边界条件的直接转移。例如，点 P(k, j)为非正则内点，其最接近的边界点为 Q点，则有上式可以看作是用

9、零次插值得到非正则内点处 u的近似值，容易求出，其截断误差为 O(h )。 将上式代入，方程个数即与未知数 个数相等。(2)线性插值这种方案是通过用同一条网格线上与点P相邻的边界点与内点作线性插值得到非正则内点 P(k, j)处U值的近似。由点R与T的线性插值确定U(P)的近似值U k, j ,得2其中d RP ,其截断误差为O(h )。将其与方程相近似的差分程联立，得到方程个数与未知数个数相等的方程组，求解此方程组可得Poisson方程第一边值问题的数值解。上面所给出的差分格式称为五点菱形格式，实际计算时经常取h ,此时五点菱形格式可化为简记为 uk,jfk,j其中 uk,j uk 1,j

10、uk 1,j例1用五点菱形格式求解拉普拉斯(其中(x, y)0 x, yuk,j 1uk,j 14uk,j0Laplace )方程第一边值问题_1(0,0) (1,0)(2,0)(3,0)解网格中有四个内点，均为正则内点。由五点菱形格式，得方 程组h2(U2,1U0,1U1 ,2U1, 04u1,1 )5 (U31U1,1 U2,2 U2,04U2,1 )0h2(U3,2U0,2 U1，3 U1,14U1,2)0 代入边界条件U1,2 U2,3 U2,14u2,2 )0U1,0U0,1U1,3U3,1其解为,16 ig I, u 9,10 ig 一, u 9,25|gm,u,37ig 一, u

11、 92,00,22,33,2U1,10.2756919 ,ig25,13喝,341g万,40 lg3U2,10.4603488u1,20.3467842, u2,20.5080467当h 时，对1卜2 (Uk 1,jUk 1,j Uk,j 1 Uk,j 1 4uk,j) fk,j禾用点(k,j), (k 1, j 1), (k 1, j 1)构造的差分格式，称为五点矩形格式，简记为. 2 口 Uk,jfk,jh其中口uk,j uk 1,j 1uk 1,j 1uk 1,j 1uk 1,j 14uk,j ,其截断误差为2五点菱形格式与矩形格式的截断误差均为O(h ),称它们具有二阶精度。如果用更多

12、的点构造差分格式，其截断误差的阶数可 以提高，如利用菱形格式及矩形格式所涉及的所有节点构造出的 九点格式就是具有四阶精度的差分格式。 3抛物型方程的差分解法以一维热传导方程为基本模型讨论适用于抛物型方程定解问题的几种差分格式。3.1差分格式的建立首先对xt平面进行网格剖分。分别取h,为x方向与t方向的步长，用两族平行直线x xkkh(k 0, 1, 2,L ),ttjj (j 0,1,2 )，将 xt 平面剖分成矩形网格，节点为(xk,tj)(k 0, 1, 2,L ,j 0,1,2L )。为简便，记(k, j) (xk,tj) , u(k, j) u(xk,tj),k (xk), g1jgt

13、j),g2jg2 (tj ),1 j1 (tj ),2j 2(tj ) o(一)微分方程的差分近似在网格内点(k, j)处，对 一J分别采用向前、向后及中心差商公式一维热传导方程可分别表示为由此得到一维热传导方程的不同差分近似2、上述差分方程所用到的节点各不相同。其截断误差分别为O( h ),222、0( h )和0( h )。因此，它们都与一维热传导方程相容。 如果将式1()中的u k, j用2 (Uk 1Uk，J 1 )代替，则可得到又一种差分近似差分方程用到四个节点。由Taylor公式容易得出,2. 2 x故其的截断误差为0( h )h 2。因而不是对任意的h,0,此差分方程都能逼近热传

14、导方程(a 0)仅当 o(h)时，才成立。综上可知，用不同的差商公式可以得到微分方程的不同的差 分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳 定性。前面三个方程都具有相容性，而此方程则要在一定条件下 才具有相容性。(二)初、边值条件的处理为用差分方法求解定解问题初值问题初边值问题还需对定解条件进行离散化。对初始条件及第一类边界条件，可直接得到 TOC o 1-5 h z lT苴中n, m 一八十h对第二、三类边界条件需用差分近似。下面介绍两种较简单的处理方法。u(1)在左边界(x 0 )处用向前差商近似偏导数在右边界(xl )处用向后差商近似则得边界条件的差分近似为其截断误差为 0

15、(h)u(2)用中心差商近似，即x则得边界条件的差分近似为2其截断误差为o(h )0误差的阶数提高了，但出现定解区域外的节点(1, j)和(n 1, j),这就需要将解拓展到定解区域外。可以通过用内节点上的U值插值求出u 1,j和un 1,j ,也可以假定热传导方程在边界上也成立，将差分方程扩展到边界节点上，由此消去 u i,j和 u n 1, j。(三)几种常用的差分格式以热传导方程的初边值问题 为例给出几种常用的差分格式。(1)古典显式格式a令r . 2 ,贝Uuk,j 1 uk,j uk 1, j 2uk, j uk 1, ja2 0h2可改写成uk,j 1 ru k 1,j(1 2r)

16、uk,j ru k 1,j将其与初始条件及第一类边界条件结合，我们得到求解此问题的一种差分格式由于第0层(j0)上节点处的u值已知(uk,0k),由此即可算出u在第一层(j1)上节点处的近似值U k,1。重复使用此式，可以逐层计算出所有的Uk,j ,因此此差分格式称为古典显式格式。又因式中只出现相邻两个时间层的节点，故此式是二层显式格式。(2)古典隐式格式将式Uk,j Uk,j 1Uk 1,j a2Uk, jUk 1,jh2整理并与初始条件及第一类边界条件式联立，得差分格式如下a其中r 卜2 。虽然第0层上的U值仍为已知，但不能由上式直接计算以上各层节点上的值 Uk,j ,必须通过解下列线性方

17、程组才能由Uk,j计算Uk,j 1 ,故此差分格式称为古典隐式格式。此方程组是三对角方程组，且系数矩阵严格对角占优，故解存在唯一。(3) Richardson 格式Richardson格式是将式整理后与初始条件及第一类边界条件式联立。其计算公式为这种差分格式中所涉及的节点出现在j 1, j, j 1三层上，故为三层显式格式。Richardson格式是一种完全不稳定的差分格式，因此它在实际 计算中是不能采用的。(4)杜福特-弗兰克尔(DoFort-Frankel )格式DoFort-Frankel 格式也是三层显式格式，它是由式与初始条件及第一类边界条件式结合得到的。具体形式如下：用这种格式求解

18、时，除了第0层上的值U k,0由初值条件得到，必须先用二层格式求出第1层上的值Uk,1 ,然后再按上式逐层计算 Uk,j(j 2,3, ,m)。(5)六点隐式格式对二阶中心差商公式2在点(k,j 1)与点(k, j) 处的二阶中心差商的平如果用 Y 2 xu均值近似2u2x(i 1)(k，J 2)处的值，即(k i )同时在点(k, j 2)处的值也用中心差商近似，即这样又得到热传导方程的一种差分近似22、其截断误差为O(h ),将上式与初始条件及第一类边界条件式联立并整理，得差分格式此格式涉及到六个节点，它又是隐式格式，故称为六点隐式格式。与古典隐式格式类似，用六点格式由第 j层的值u k,

19、 j计算第j 1层的值uk,j 1时，需求解三对角方程组此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。例2用古典显式格式求初边值问题 TOC o 1-5 h z 的数值解，取h 1,0.5 0a 1,r ar0.5(x)x2解这里 h h2，(x)x ,g) 0,g2(t) 9。由格式可得到将初值u k ,0代入上式，即可算出将边界条件u0,10,u3,19及上述结果代入又可求得如此逐层计算，得全部节点上的数值解为 4双曲型方程的差分解法对二阶波动方程_u_u如果令v1f , v2x ,则方程可化成一阶线性双曲型方程组记v (Vi, V2 )T ,则方程组可表成矩阵形式矩阵 A有两个不同的特征值a,故存在非奇异矩阵 P,使得作变换w Pv (wi,w2),方程组可化为w wtx方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。4.1