-偏微分方程数值解法

1、精品第十章偏微分方程数值解法偏微分方程问题，其求解十分困难。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。§1差分方法的基本概念1.1几类偏微分方程的定解问题椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程2uU 2 x2u2yf (x,y)特别地，当 f （x, y） 0时，即为拉普拉斯（Laplace ）方程，又称为调和方程22u u cU 220x yPoisson方程的第一边值问题为22(x, y)u u1 1 f(x, y) x yu(x, y) (x,y)(x, y)-可编

2、辑-精品其中 为以 为边界的有界区域，为分段光滑曲线，称为定解区域，f (x, y), (x, y)分别为 ，上的已知连续函数。第二类和第三类边界条件可统一表示为u/ u (x,y)(x,y)n其中n为边界的外法线方向。当0时为第二类边界条件，0时为第三类边界条件。抛物型方程：其最简单的形式为一维热传导方程2uua2 0 (a 0)tx方程可以有两种不同类型的定解问题：初值问题2u ua 0 t 0, x t xu(x,0)(x)x初边值问题-可编辑-精品2-u aur 0 0 t T, 0 x ltx2u(x,0)(x)0 xlu(0,t)g1(t),u(l,t)g2(t)0 tT其中 (x

3、), g1(t), g2(t)为已知函数，且满足连接条件(0)gi(0),(i) g2(0)边界条件 u(0,t)gi(t),u(l,t)g2(t)称为第一类边界条件。第二类和第三类边界条件为i(t)u x 0 gi(t)2(t)u x l g 2 (t)其中 i(t) 0, 2(t) 0。当 l(t)2 (t)0时，为第二类边界条件,-可编辑-精品否则称为第三类边界条件双曲型方程：最简单形式为一阶双曲型方程uuC a0tx物理中常见的一维振动与波动问题可用二阶波动方程2u2x描述，它是双曲型方程的典型形式。方程的初值问题为22t 0,u 2 u.2 a 2t Xu(x,0)(x)(x)边界条

4、件一般也有三类，最简单的初边值问题为-可编辑-精品u 2 u2 a 2 00 t T, 0 x l t xu(x,0)(x),U(x) 0 x lt t 0u(0,t) gi(t), u(l,t) g2(t)0 t T1.2差分方法的基本概念差分方法又称为有限差分方法或网格法，是求偏微分方程定 解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连 续变化区域用有限离散点(网格点)集代替；将问题中出现的连 续变量的函数用定义在网格点上离散变量的函数代替；通过用网 格点上函数的差商代替导数，将含连续变量的偏微分方程定解问 题化成只含有限个未知数的代数方程组(称

5、为差分格式)。如果 差分格式有解，且当网格无限变小时其解收敛于原微分方程定解 问题的解，则差分格式的解就作为原问题的近似解(数值解)。 因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：(1)选取网格；(2)对微分方程及定解条件选择差分近似，列出差分格式；(3)求解差分格式；(4)讨论差分格式解对于微分方程解的收敛性及误差估计。-可编辑-精品下面，用一个简单的例子来说明用差分方法求解偏微分方程问题的一般过程及差分方法的基本概念。设有一阶双曲型方程初值问题。t 0,u -a 0xu(x,0)(x)(1)选取网格:a-,-2h-h0 h2r3hD ()(,t)x,t 0作网检；剖首先对定解

6、区域分，最简单常用一种网格是用两族分别平行于x轴与轴的等距直线x xkkht tjj (k 0, 1, 2L ,j0,1,2,L )将D分成许多小矩形区域。这些直线称为网格线，其交点称为网格点，也称为节点,h和 分别称作 x 方向和t方向的步长。这种网格称为矩形网格。-可编辑-精品uu(4 i,tjxtj)uu(Xk,tj 1)u(Xk,tj)h ,ha%")u(Xk,tj),(2)对微分方程及定解条件选择差分近似，列出差分格式:如果用向前差商表示一阶偏导数，即t (Xk,tj)其中01, 21。u、灯a方程t在下点(X k , t j )处口表小为u(Xk,tj i) u(Xk,t

7、j)_2Ut2(Xk,tjuc0Xu(Xk i,tj) u(Xk,tj)"ut2 (Xk,tjR(Xk,tj)(其中 u(Xk,0)a h2 )”2乂ih,tj)k 0, 1, 2, ,j 0,1,2,)(Xk)(k0,1, 2,L)。由于兰-可编辑-精品h, 足够小时，在式中略去 R(xk,tj),j,就得到一个与方程相近似的差分方程Uk 1,jah此处，Uk,j可看作是问题的解在节点(xkj )处的近似值。同初值条件Uk,0(Xk)(k 0, 1, 2,)结合，就得到求问题的数值解的差分格式。R(Xk,tj) 2ut2(Xk,tj2 ) ux2(Xkih,tj)O( h)称为差分

8、方程的截断误差。如果一个差分方程的截断误差为R O( q hp),则称差分方 程对t是q阶精度，对x是p阶精度的。显然，截断误差的阶数 越大，差分方程对微分方程的逼近越好。若网格步长趋于0时，差分方程的截断误差也趋于0,则称-可编辑-精品差分方程与相应的微分方程是相容的。这是用差分方法求解偏微分方程问题的必要条件。如果当网格步长趋于0时，差分格式的解收敛到相应微分方程定解问题的解，则称这种差分格式是收敛的。§2椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法。2.1差分格式的建立(x,y)考虑Poisson方程的第一边值问题f(x, y)-

9、可编辑-精品R (Xk,yj)Xk kh,yj, k, j为整数。定解区域内部的节点称为内点，记内点集 R为 h 。边界 与网格线的交点称为边界点，边界点全体记为。与节点（xk , y j ）沿X方向或y方向只差一个步长的点（xk i,yj）和（xk,yj i）称为节点(Xk,yj)的相邻节点。如果一个内点的四个相邻节点均属于则内点，正内点的全体记为,至少有一个相邻节点不属于的内点称为非正则内点，非正则内点的全体记为。问题是要求出第一边值问题在全体内点上的数值解。为简便，记(k,j) (Xk,yj) , u(k,j) u(xk,yj),fk,j f (xk,yj)0 对正则内点（k, j），、

10、，,由二阶中心差冏公式-可编辑-精品2 u2X2 u2y(k,j)(k,j)Poisson 方程(k,j)处可表示为u(k 1,j) u(k, j)u(k,j) u(k 1,j)u(k 1,j) 2u(k, j) u(k 1,j)h2u(k,j 1) 2u(k,j) u(k,j 1)2u2X2u2yu(k 1,j) 2u(k,j)u(k 1,j)其中h2h2R(k,j) 12uX4)(Xkh,yj)£12h2u(t12 xu：? (Xkf(x，y)在点u(k,j 1) 2u(k, j) u(k,212uX4)(xk,yj为其截断误差表示式，略去 R(k,j), 即得与方程相近似的差分

11、方程uk 1,j 2uk, j uk 1,jh2uk,j 1 2uk,j uk,j 12-可编辑-2 ) O(h2fk,j精品式中方程的个数等于正则内点的个数，而未知数 uk,J则除了包含 正则内点处解u的近似值外，还包含一些非正则内点处 u的近 似值，因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按上式给出，需要利用边界条件得到。边界条件的处理可以有各种方案，下面介绍较简单的两种。(1)直接转移用最接近非正则内点的边界点上的u值作为该点上u值的 近似，这就是边界条件的直接转移。例如，点P(k, j)为非正则内点，其最接近的边界点为 Q 点，则有uk,ju(Q) (

12、Q)(k,j)上式可以看作是用零次插值得到非正则内点处 u的近似值，容易 求出，其截断误差为 O(h ) 0将上式代入，方程个数即与未知数 个数相等。(2)线性插值这种方案是通过用同一条网格线上与点 P相邻的边界点与 内点作线性插值得到非正则内点 P(k, j)处u值的近似。由点R与T的线性插值确定 u(P) 的近似值u k , j ,得-可编辑-精品Uk,jh-d(R)(T)其中dC-八 2 RP，其截断误差为O(h )。将其与方程相近似的差分程联立，得到方程个数与未知数个数相等的方程组，求解此方程组可得Poisson方程第一边值问题的数值解。上面所给出的差分格式称为五点菱形格式,Uk,j

13、1fk,jUk 1,j 2uk, j Uk 1, j Uk,j 1 2uk, jh2r实际计算时经常取h,此时五点菱形格式可化为-1(.卜2 (Uk 1,j Uk 1,juk,j1 Uk,j14uk,j )fk,j简记为Uk 1,jUk,jUk,j 1fk,jUk,j 14uk,j0例1用五点菱形格式求解拉普拉斯(Laplace )方程第一边值问题-可编辑-精品其中f (x, y)(x, y)u(x, y)ig(1 x)2y2( x, y)0 x, y（0,0） （10)(2,0)(3,0)解网格中有四个内点，均为正则内点C由五点菱形为式，得方程组记伍1u0,1Ui,2U1,04ui,i)土伍

14、1U1,1U2,2U2,04u2,1 )1h21h2(U2,2(U3,2U0,2U1,2U1,3U2,3U1,1U2,1-可编辑-4ui,2 )4u2,2 )0 代入边界条件精品u1,0U0,1U1,3U3,1161g3, u,101gj,u,25ig I, u 9,371g J,u2,00,22,33,2其解为U1,10.275691925 ig 9,131g 8igigu1,20.3467842, u2,234940u2,10.46034880.5080467Li卜2 (Uk 1,jUk 1,jUk,jUk,j 14uk,j)fk,j利用点(k, j), (k 1, j1), (k 1,

15、j 1)构造的差分格式，称为五点矩形格式，简记为uk,jfk,juk,j uk 1,j 1Uk 1, j 1 Uk 1,j 1Uk 1,j 1 4uk,j ,其截断误差为-可编辑-精品2444h uu u _R(k,j)4 6244 O(h12 xx y y *,j)2、五点菱形格式与矩形格式的截断误差均为 0(h),称它们具有二阶精度。如果用更多的点构造差分格式，其截断误差的阶数可以提高，如利用菱形格式及矩形格式所涉及的所有节点构造出的九点格式就是具有四阶精度的差分格式。§3抛物型方程的差分解法以一维热传导方程2u ua10 (a 0)t x为基本模型讨论适用于抛物型方程定解问题的

16、几种差分格式。3.1差分格式的建立首先对 xt平面进行网格剖分。分别取h,为X方向与t方向的步长，用两族平行直线x xk kh(k 0, 1, 2，L )t tj j (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)(xk)gj g1(tj)g2j g2(tj), 1j1 (tj ), 2j 2(tj )-可编辑-精品(一)微分方程的差分近似在网格内点(k, j)处，对u 分别采用向前、向后及中心差商公式u(k,j 1) u(k,j)O()、(k

17、,j)ut (k,j)ut (k,j)u(k,j) u(k,j 1)O()u(k,j 1) u(k,j 1)2O( 2)一维热传导方程2u ca2-0x(a 0)可分别表示为u(k,j 1) u(k, j)au(k 1,j)2u(k,j) u(k 1,j) O( h2u(k,j 1) u(k,j 1) au(k 1,j) 2u(k,j) u(k 1,j)h2O( Iu(k,j 1) u(k,j 1) _u(k 1,j) 2u(k,j) u(k 1,j)ah2O( 2 h：由此得到一维热传导方程的不同差分近似-可编辑-精品uk,j 1 uk, juk 1,j a2uk, j uk 1,j h1u

18、k,j uk,j 1uk 1,j a2uk, j uk 1,juk,j 1 uk,j 12uk 1, j a2uk, j uk 1, j h上述差分方程所用到的节点各不相同其截断误差分别为0( h2),0( h2)和0( 2 h2)因此，它们都与一维热传导方程相容。如果将式u(k,j 1) u(k,j 1)u(k 1,j) 2u(k,j) u(k 1,j)h20(21()中的u k, j用2 (uk,J 1 uk" 1)代替，则可得到又一种差分近似uk,j 1uk,j 1uk 1, j uk,j 1 uk,j 1uk 1,ja2 0h2差分方程用到四个节点。由Taylor公式容易得出

19、-可编辑-精品12、uk,j 2(Uk,j 1 uk,j 1 ) O()故其的截断误差为0(h2)h 2。因而不是对任意的h,0,此差分方程都能逼近热传导方程2- a 0-0 (a 0)t x2仅当 o(h)时，才成立。综上可知，用不同的差商公式可以得到微分方程的不同的差 分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳 定性。前面三个方程都具有相容性，而此方程则要在一定条件下 才具有相容性。(二)初、边值条件的处理为用差分方法求解定解问题初值问题2u ua r 0 t 0, x tx2u(x,0)(x)x初边值问题-可编辑-精品2u2Xt T, 0 x lu(x,0)(x)u(0,

20、t)gi(t),u(l,t)g2(t)还需对定解条件进行离散化。对初始条件及第一类边界条件,可直接得到u k,0u(Xk，0)(k0,1,L 或 k0,1,L ,n)0, ju(0,tj)giju(l,tjg2j0,1, m1)其中n对第二、三类边界条件(t)ug i (t)2 (t )ug 2 (t)需用差分近似。下面介绍两种较简单的处理方法。-可编辑-精品(1)在左边界(xu0 )处用向前差商近似偏导数二 ,x在右边界(xl )处用向后差商近似(0, j)u(1, j) u(0, j) O(h)hu(n, j) u (n 1, j)O (h)(j 0,1,m)则得边界条件的差分近似为u1,

21、j u0,jhu0,j g1ju n,j u n 1, jh2 ju n, j 9 2 j(j 0,1,m)其截断误差为 O(h)(2)用中心差商近似u(1,j)(0, j)u(2h5 O(h2)u(n 1,j) u(n 1,j)(n,j)2hO(h2)-可编辑-精品(j 0,1,m)则得边界条件的差分近似为ui,j u i,j2hi ju 0, j g i j2 j Un, j g 2 j(j 0,1,m)2、其截断误差为O(h ) 0误差的阶数提高了，但出现定解区域外的节点(1, j)和(n 1, j),这就需要将解拓展到定解区域外。可以通过用内节点上的U值插值求出U 1,j和Un 1,j

22、 ,也可以假定热传导方程在边界上也成立，将差分方程扩展到边界节点上，由此消去u 1,j和 U n 1, j。(三)几种常用的差分格式2Ua2xU(x,0)(x)u(0, t) g4t),以热传导方程的初边值问题0 0 t T, 0 x l0 x lU(l,t)g2(t)0 t T为例给出几种常用的差分格式(1)古典显式格式-可编辑-精品uk,j 1 u k,j uk 1, j 2uk, j uk 1, ja 0h2可改写成Uk,j i ru k i,j (1 2r)Uk,j ru k i,j将其与初始条件及第一类边界条件Uk,0 U(Xk，0) k(k 0, 1,L 或 k 0,1,L ,n)

23、U0,ju(0,tj )g1jUn,ju(1 ,tj )g2j(j 0,1, m 1)结合，我们得到求解此问题的一种差分格式Uk,j 1rUk 1,j(12r)Uk,jrUk 1(k 1,2, ,n 1,j 01, ,mUk,0 kU0, jg1j , Un, j g2j(k 0,1, ,n)(j 0,1, ,m)由于第U在0层(j 0)上节点处的U值已知(uk,0k) ,由此即可算出-可编辑-精品第一层(j1)上节点处的近似值Uk,1 o重复使用此式，可以逐层计算出所有的Uk,j ,因此此差分格式称为古典显式格式。又因式中只出现相邻两个时间层的节点，故此式是二层显式格式。(2)古典隐式格式将

24、式uk,jUk,j 1 Uk 1,j 2uk,juk 1,ja； 0 h2整理并与初始条件及第一类边界条件式联立，得差分格式如下Uk,j 1 -j r(Uk 1,j 1 2%j1 Uk 1,j 1) (k 1,2L ,n 1,j 0Uk,。 k(k 0,1L,n)% *, Un,jg2j(j 0,1L ,m)a其中r 卜2 。虽然第0层上的U 值仍为已知，但不能由上式直接计算以上各层节点上的值 Uk,j ,必须通过解下列线性方程组Uk,j1r(Uk1,j 12Uk,j1Uk1,j 1)rUk1,j 1(1 2r)Uk,j 1- 1,j-可编辑-精品1 2r2ru1,j 1u1,j rg2r2r

25、0,1,才能由uk,j计算uk,jU2, j 1Mu2, jMun2, j 1i,j i1)1 ,故此差分格式称为古典隐式格式。此方程组是三对角方程组，且系数矩阵严格对角占优，故解存在唯(3) Richardson 格式Richardson 格式是将式u k, j 1 u k, j 1k 1, j 2uk, j u k 1, jh"整理后与初始条件及第一类边界条件式联立其计算公式为un 2,un 1,jr(uk,j 1uk,j 12r(uk 1,j2uk,juk 1,j )(k1,2, ,n 1,uk,0(k0,1,n)u0,j91j, un,j g2j(j0,1,m)这种差分格式中

26、所涉及的节点出现在j1, j, j1三层上，故为三层显式格式。Richardson 格式是一种完全不稳定的差分格式，因此它在实际-可编辑-精品uk,j 1 uk 1, jh22rUk,j1 2r(uk 1,j uk 1,j )uk,j1 2r(k12,nuk,0(k01,n)U0,jg1j,un,j g2j(j01,m)计算中是不能采用的。(4)杜福特-弗兰克尔(DoFort-Frankel )格式DoFort-Frankel格式也是三层显式格式，它是由式uk,j 1uk, j 1uk 1, j uk, j 1a与初始条件及第一类边界条件式结合得到的。具体形式如下:用这种格式求解时，除了第0层

27、上的值u k,0由初值条件得到，必须先用二层格式求出第1层上的值uk,1 ,然后再按上式逐层计算 uk,j(j 2,3,m)o对二阶中心差商公式2u2xu(k 1,j) 2u(k,j)u(k 1,j)(k, j)h2-可编辑-精品如果用2u2x在点(k, j 1)与点(k, j) 处的二阶中心差商的平均值(k, j2)处的值,即2u近似一彳在x(k,juk 1,j 1 2uk,j 1, uk 1,j 1 uk 1,j 2uk, j u 2h2同时ut 在点(k,j1) A / R2,处的值也用中心差冏近似，即(k,j 2)这样又得到热传导方程的一种差分近似uk,j 1uk,ja /2 (u k

28、 1, j 1 2uk, j 1 , u k 1, j 1 uk 1, j 2u k, j 2h22、其截断误差为O(h ),将上式与初始条件及第一类边界条件式联立并整理，得差分格式-可编辑-精品uk,j 1ruk, j 2 (uk 1,j 12uk,j 1uk 1,j、r /1)2 (uk 1，j 2uk,j(k1,2,L ,n1)uk,0(k0,1L ,n)u0,jg1j, un,jg2j(j 0,1L ,m)此格式涉及到六个节点，它又是隐式格式,故称为六点隐式格式。与古典隐式格式类似，用六点格式由第j层的值u计算第层的值uk,j 1时，需求解三对角方程组1 r r/2. r/2 1 rr

29、/2u1,j 1u2,j 1Mun 2,j 1r/21 r r/2Un 1,j 1(j 0,1, m1)此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。例2用古典显式格式求初边值问题-可编辑-r/2 1 r精品3,u(x,0)u(0,t)0,u(3,t)的数值解，取h1,0.5。“a解 这里1,ra2 0.5 h2(x)gi(t)0,g2(t) 9。由格式uk,ji,j(i2r)uk,jruki,j(k1,2,n 1,可得到uk,juk,0u0,j将初值uk,0u0,j(k0,i,n)gij,un, jg2j(j0,i,m)0.5(uk2xkk20, u3,jk ,0代入上式,ui,i0.

30、5i,juk i,j )即可算出(u2,0u0,0)0.5-可编辑-(k(k(j(4i,2, j0,i,5)0,i,2,3)0,i, ,6)0) 2精品u210.5 (u30 u10) 0.5 (9 1) 5将边界条件U0,10, u319及上述结果代入又可求得U0,20ul ,225U2,25.5, U329如此逐层计算，得全部节点上的数值解为U0,30, U1,3 2.75, U2,3 5.75, u、U0,40, L L§4双曲型方程的差分解法对二阶波动方程如果令V1V2则方程可化成一阶线性双曲型方程V1V2V1v(V1, V2 )T,则方程组可表成矩阵形式-可编辑-精品v0 a2VA vA矩阵A有两个不同的特征值PAPt10x xa,故存在非奇异矩阵 P,使得a 00 a作变换W Pv (W|,W2),方程组可化为w wtx方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。4.1 几种简单的差分格式考虑一阶双曲型方程的初值问题uua 0 t 0, x txu(x,0)(x)x将xt平面剖分成矩形网格，取x方向步长为h, t方向步长为网格线