偏微分方程的数值解法

1、第十六章 偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial differential equation,简称PDE 。偏微分方程问题,其求解是十分困难的。除少数特殊情况外,绝大多数情况均难以求出精确解。因此,近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下:(,x x x y y y x y A B C f x y += (16.1.1 式中,A B C 是常数,称为拟线性(quasilinear数。通常,存在3种拟线性方

2、程: 双曲型(hyperbolic方程:240B AC -; 抛物线型(parabolic方程:240B AC -=; 椭圆型(ellliptic方程:240B AC -+=-+ (16.1.4边界条件一般有三类,最简单的初边值问题为:2222212000,0(,0(0,(,(,(0(t u ua t T x l t x u x lu t g t u l t g t t T ux x t =- (16.1.8 方程可以有两种不同类型的定解问题:(1 初值问题:2200,(,0(u ua t x t xu x x x -=-+=-+(16.1.6(2 初边值问题:221200,0(,0(0(0,

3、(,(,(0u ua t T x l t x u x x x l u t g t u l t g t t T-=-+= (16.2.1选取网格: 图16.2.1 差分示意图首先对定解区域(,0D x t x t =-+作网格剖分,最简单常用一种网格是用两族分别平行于x 轴与t 轴的等距直线:k x x k h =,(0,1,2,0,1,2,j t t j k j = (16.2.2将D 分成许多小矩形区域。这些直线称为网格线,其交点称为网格点,也称为节点,h 和分别称作x 方向和t 方向的步长。这种网格称为矩形网格。(1 对微分方程及定解条件选择差分近似,列出差分格式。如果用向前差商表示一阶偏

4、导数,即:2211(,12(,(,(,(,2(,(,(,2k j k j k j k j k j x x t k j k j k j t x t u x t u x t u h u x h t x h u x t u x t u u x t t +-=-+-=-+ (16.2.3其中120,1-+=-+=-+ 其中:101(0200x x x x ,用差分格式求其数值解,(1,2,3,4k j u j =,取12r h=。 解:记(0,1,2,k x kh k =,由初始条件:11,2,1(0201,2,k k k x k k =-= 按差分格式:,1,1,0(0,1,2,0,1,2,k j

5、k j k j k j k ku u ar u u k j u +=-=计算公式为:,1,1,3122k j k j k j u u u +=-。计算结果略。如果用差分格式:,1,1,0(0,1,2,0,1,2,k j k j k j k j k ku u ar u u k j u +-=-=求解,计算公式为:,1,1,1(2k j k j k j u u u +-=+计算结果略。与准确解:11(,200x t u x t x t x 比较知,按前一个差分格式所求得的数值解不收敛到初值问题的解,而后一个差分格式的解收敛到原问题的解。16.3.2 波动方程的差分格式对二阶波动方程:22222u

6、u a t x = (16.3.1 如果令1u v t =,2uv x=,则方程可化成一阶线性双曲型方程组: 21221v v a t x v v tx = (16.3.2 记T 12v (,v v =,则方程组可表成矩阵形式:20v v v10a A t x x = (16.3.3 矩阵A 有两个不同的特征值a =,故存在非奇异矩阵P ,使得:100a PAP a -=- (16.3.4 作变换12w v (,T P w w =,方程组可化为:w wt x= (16.3.5 方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。下面给出如下波动方程和边界条件的差

7、分格式: 2(,(,0,0t t x x u x t c u xy xa t b= (16.3.6 (0,0,(,00(,0(0(,0(0tu t u a t t bu x f t x a u x g x x a= (16.4.1 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。16.4.1 差分格式的建立首先对x t -平面进行网格剖分。分别取,h 为x 方向与t 方向的步长,用两族平行直线(0,1,2,k x x kh k =,(0,1,2j t t j j =,将x t -平面剖分成矩形网格,节点为:(,(0,1,2,0,1,2k j x t k j =。为简便,记:(,(,k j

8、 k j x t =,(,(,k j u k j u x t =,(k k x =,(k k x =,22(j j g g t =,11(j j t =,22(j j t =。16.4.2 微分方程的差分近似在网格内点(,k j 处,对u t分别采用向前、向后及中心差商公式: (,(,2(,(,1(,(,(,1(,1(,1(2kj kj k j uu k j u k j O t uu k j u k j O t uu k j u k j O t +-=+-=+-=+ (16.4.2 一维热传导方程可分别表示为:22(,1(,(1,2(,(1,(u k j u k j u k j u k j u

9、 k j a O h h +-+-+-=+ (16.4.3 22(,1(,1(1,2(,(1,(u k j u k j u k j u k j u k j a O h h+-+-+-=+ (16.4.4 222(,1(,1(1,2(,(1,(2u k j u k j u k j u k j u k j a O h h +-+-+-=+ (16.4.5 由此得到一维热传导方程的不同差分近似:,1,1,1,220k j k j k j k j k j u u u u u a h +-+-= (16.4.6 ,11,1,220k j k j k j k j k ju u u u u a h -+-+

10、-= (16.4.7,1,11,1,220k j k j k j k j k ju u u u u a h +-+-+-= (16.4.8上述差分方程所用到的节点各不相同。其截断误差分别为2(O h +,2(O h +和22(O h +。因此,它们都与一维热传导方程相容。 如果将式:222(,1(,1(1,2(,(1,(2u k j u k j u k j u k j u k j a O h h+-+-+-=+中的,k j u 用,1,11(2k j k j u u +-+代替,则可得到又一种差分近似: ,1,11,1,11,202k j k jk j k j k j k j u u u uu

11、 u a h +-+-+-= (16.4.9差分方程用到四个节点。由Taylor 公式容易得出:2,1,11(2k j k j k j u u u O +-=+ (16.4.10 故其的截断误差为2222(O h O h + 。因而不是对任意的,0h ,此差分方程都能逼近热传导方程:220,(0u u a a t x-=,仅当(o h =时,才成立。 综上可知,用不同的差商公式可以得到微分方程的不同的差分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳定性。前面三个方程都具有相容性,而此方程则要在一定条件下才具有相容性。16.4.3 初、边值条件的处理为用差分方法求解定解问题初值问题:2

12、200,(,0(u u a t x t xu x x x -=-+=-+ (16.4.11初边值问题:221200,0(,0(0(0,(,(,(0u u a t T x l t x u x x x l u t g t u l t g t t T -=(16.4.12 还需对定解条件进行离散化。对初始条件及第一类边界条件,可直接得到:,0(,0k k k u u x =,(0,1, 0,1,k k n =或 0,1,2(0,(,j j j n j j j u u t g u u l t g =,(0,1,1j m =- 式中:,l T n m h =。 对第二、三类边界条件:10122(x x

13、l u t u g t x u t u g t x =-=+= 0t T (16.4.13 需用差分近似。下面介绍两种较简单的处理方法。在左边界(0x =处用向前差商近似偏导数u x ,在右边界(x l =处用向后差商近似u x,即: (0,(,(1,(0,(,(1,(j n j u u j u j O h x h u u n j u n j O h x h-=+-=+,(0,1,j m = (16.4.14 则得边界条件的差分近似为:1,0,10,1,1,2,2j j j j j n j n j j n j j u u u g h u u u g h -=-+=,(0,1,j m = (16

14、.4.15其截断误差为(O h 。(2 用中心差商近似u x,即: 2(0,2(,(1,(1,(2(1,(1,(2j n j u u j u j O h x hu u n j u n j O h x h -=+-=+,(0,1,j m = (16.4.16则得边界条件的差分近似为:1,1,10,11,1,2,22j j j j j n j n j j n j j u u u g h u u u g h -+-=-+=, (0,1,j m = (16.4.17其截断误差为2(O h 。误差的阶数提高了,但出现定解区域外的节点(1,j -和(1,n j +,这就需要将解拓展到定解区域外。可以通过用

15、内节点上的u 值插值求出1,j u -和1,n j u +,也可以假定热传导方程在边界上也成立,将差分方程扩展到边界节点上,由此消去j u ,1-和j n u ,1+。16.4.4 几种常用的差分格式以热传导方程的初边值问题:221200,0(,0(0(0,(,(,(0u u a t T x l t x u x x x l u t g t u l t g t t T -=(16.4.18 为例给出几种常用的差分格式。(1 古典显式格式 令2a r h =,则: ,1,1,1,220k j k j k j k j k j u u u u u a h +-+-= (16.4.19可改写成:,11,

16、1,(12k j k j k j k j u ru r u ru +-=+-+ 。将其与初始条件及第一类边界条件: ,0(,0k k k u u x =,(0,1, 0,1,k k n =或0,1,2(0,(,j j j n j j ju u t gu u l t g =,(0,1,1j m =-结合,我们得到求解此问题的一种差分格式:,11,1,00,1,2(12(1,2,1,0,1,1(0,1,(0,1,k j k j k j k j k kjj n j j u ru r u ru k n j m u k n u g u gj m +-=+-+=-=-= (16.4.20由于第0层(0j

17、=上节点处的u 值已知(0,k k u =,由此即可算出u 在第一层(1j =上节点处的近似值,1k u 。重复使用此式,可以逐层计算出所有的,k j u ,因此此差分格式称为典显式格式。又因式中只出现相邻两个时间层的节点,故此式是二层显式格式。(2 古典隐式格式 将式:,11,1,220k j k j k j k j k ju u u u u a h-+-+-= (16.4.21 整理并与初始条件及第一类边界条件式联立,得差分格式如下:,1,1,1,11,1,00,1,2(2(1,2,1,0,1,1(0,1,(0,1,k j k j k j kj k j k kjj n j j u u r

18、u u u k n j m u k n u g u gj m +-+=+-+=-=-= (16.4.22其中:2a r h=。虽然第0层上的u 值仍为已知,但不能由上式直接计算以上各层节点上的值,k j u ,必须通过解下列线性方程组:1,11,112,12,2,12,1,11,2112012001212j j j j j n j n j n j n j j r r u u rg r r ru u u u u u rg rr r r r +-+-+-+-+-+-=+-+-+(16.4.23 式中,(0,1,1j m =-。 才能由,k j u 计算,1k j u +,故此差分格式称为古典隐式格

19、式。此方程组是三对角方程组,且系数矩阵严格对角占优,故解存在唯一。(3 Richardson 格式 Richardson 格式是将式:,1,11,1,220k j k jk jk j k ju u uu u ah +-+-+-=整理后与初始条件及第一类边界条件式联立。其计算公式为:,1,11,1,00,1,22(2(1,2,1,1,2,1(0,1,(0,1,k j k j k j k j k j k k jj n j j u u r u u u k n j m u k n u g u g j m +-+-=+-+=-=-= 6.4.24这种差分格式中所涉及的节点出现在1,1+-j j j 三层

20、上,故为三层显式格式。Richardson 格式是一种完全不稳定的差分格式,因此它在实际计算中是不能采用的。(4 杜福特-弗兰克尔 (DoFort-Frankel 格式 DoFort-Frankel 格式也是三层显式格式,它是由式:,1,11,1,11,202k j k jk jk j k j k ju u uuuuah+-+-+-= 与初始条件及第一类边界条件式结合得到的。具体形式如下:,11,1,1,00,1,2212(1,2,1,1,2,11212(0,1,(0,1,k j k j k j k j k kj j n j j r r u u u u k n j m r ru k n u g

21、 u g j m +-=+=-=-+=(16.4.25用这种格式求解时,除了第0层上的值,0k u 由初值条件得到,必须先用二层格式求出第1层上的值1,k u ,然后再按上式逐层计算,(2,3,k j u j m =。(5 六点隐式格式 对二阶中心差商公式:222(,(1,2(,(1,k j u u k j u k j u k j xh +-+-=如果用22u x 在(,1k j +与(,k j 处的二阶中心差商的平均值近似22ux在1(,2k j +处的值,即:21,1,11,11,1,2221(,22,2(2k j k j k j k j k j k jk j u u u u u u uO

22、 h x h +-+-+-+-+=+同时u t 在点21,(+j k 处的值也用中心差商:(,(,1(,1(2kj u u k j u k j O t +-=+近似,即:2,1,21(,2k j k jk j u uux +-=这样又得到热传导方程的一种差分近似:,1,1,1,11,11,1,2(2,202k j k jk j k jk jk j kjk j u uau u uu u u h+-+-+-+= (16.4.26 其截断误差为(22h O +,将上式与初始条件及第一类边界条件式联立并整理,得差分格式:,1,1,1,11,11,1,00,1,2(2(222(1,2,1,0,1,1(0

23、,1,(0,1,k j k j k j k j k j k j k j k j k k j j n j jr r u u u u u u u u k n j m u k n u g u g j m +-+-=+-+-+=-=-= (14.5.27 此格式涉及到六个节点,它又是隐式格式,故称为六点隐式格式。与古典隐式格式类似,用六点格式由第j 层的值,k j u 计算第1j +层的值,1k j u +时,需求解三对角方程组:1,12,12,11,11,0,12,1,0,2,3,2,1,2,1,2,3,1/20/21/200/21/2/21(222(22(22j j n j n j j j j j

24、 j j j j j n j n j n j n j n r r u r r r u u u r r r r r r r u u u u u r u u u u r u u u u u +-+-+-+-+-+-+-+-+=+-+1,1,1,2,(222j n j n j n j n j r r u u u u +-+-+(16.4.28 式中,(0,1,1j m =-。此方程组的系数矩阵严格对角占优,故仍可用追赶法求解。例16.4.1 用古典显式格式求初边值问题。222003,03(,003(0,0,(3,903u u t x t x u x x x u t u t t -=的数值解,取1,0

25、.5h =。解:这里:21,0.5a a r h=,2(x x =,12(0,(9g t g t =。 由格式:,11,1,00,1,2(12(1,2,1,0,1,1(0,1,(0,1,k j k j k j k j k kj j n j j u ru r u ru k n j m u k n u g u gj m +-=+-+=-=-=可得到:,11,1,22,00,3,0.5(1,2,0,1,5(0,1,2,30,9(0,1,6k j k j k j k k j j u u u k j u x k k u u j +-=+=将初值,0k u 代入上式,即可算出:1,12,00,00.5(0

26、.5(402u u u =+=+= 2,13,01,00.5(0.5(915u u u =+=+=将边界条件0,10u =,3,19u =及上述结果代入又可求得:0,21,22,23,20, 2.5, 5.5,9u u u u =如此逐层计算,得全部节点上的数值解为:0,31,32,33,30,40, 2.75, 5.75,9,0,u u u u u = =16.4.5 前向差分法求解热传导方程的MATLAB 程序设(,0(u x f x =,其中,0x a ,而且,12(0,(,u t c u a t c =,其中0t b ,求解区间(,:0,0R x t x a t b =内2(,(,t

27、xx u x t c u x t =的近似解。*function U=forwdif(f,c1,c2,a,b,c,n,m%Input - f=u(x,0 as a string f % - c1=u(0,t and c2=u(a,t% - a and b right endpoints of 0,a and 0,b % - c the constant in the heat equation% - n and m number of grid points over 0,a and 0,b %Output - U solution matrix; analogous to Table 10.

28、4%Initialize parameters and Uh=a/(n-1; k=b/(m-1; r=c2*k/h2; s=1-2*r;U=zeros(n,m;%Boundary conditionsU(1,1:m=c1; U(n,1:m=c2;%Generate first rowU(2:n-1,1=feval(f,h:h:(n-2*h;%Generate remaining rows of Ufor j=2:mfor i=2:n-1U(i,j=s*U(i,j-1+r*(U(i-1,j-1+U(i+1,j-1; end endU=U;*16.4.6 Crank-Nicholson 求解热传导

29、方程的MATLAB 程序设(,0(u x f x =,其中,0x a ,而且,12(0,(,u t c u a t c =,其中0t b ,求解区间(,:0,0R x t x a t b =内2(,(,t xx u x t c u x t =的近似解。*function U=crnich(f,c1,c2,a,b,c,n,m%Input - f=u(x,0 as a string f% - c1=u(0,t and c2=u(a,t% - a and b right endpoints of 0,a and 0,b% - c the constant in the heat equation%

30、- n and m number of grid points over 0,a and 0,b %Output - U solution matrix; analogous to Table 10.5%Initialize parameters and Uh=a/(n-1;k=b/(m-1;r=c2*k/h2;s1=2+2/r;s2=2/r-2;U=zeros(n,m;%Boundary conditionsU(1,1:m=c1;U(n,1:m=c2;%Generate first rowU(2:n-1,1=feval(f,h:h:(n-2*h;%Form the diagonal and

31、off-diagonal elemnts of A and%the constant vector B and solve tridiagonal system AX=BVd(1,1:n=s1*ones(1,n;Vd(1=1;Vd(n=1;Va=-ones(1,n-1;Va(n-1=0;Vc=-ones(1,n-1;Vc(1=0;Vb(1=c1;Vb(n=c2;for j=2:mfor i=2:n-1Vb(i=U(i-1,j-1+U(i+1,j-1+s2*U(i,j-1;endX=trisys(Va,Vd,Vc,Vb; U(1:n,j=X; end U=U*16.5 椭圆型方程第一边值问题的差

32、分解法本节以Poisson 方程为基本模型讨论第一边值问题的差分方法。16.5.1 差分格式的建立考虑Poisson 方程的第一边值问题:2222(,(,(,(,(,x y u uf x y x y x yu x y x y += (16.5.1 图16.5.1取,h 分别为x 方向和y 方向的步长,如图所示,以两族平行线:k x x kh =,(,0,1,2,j y y j k j =将定解区域剖分成矩形网格。节点的全体记为: (,k j k j R x y x k h y j k j =为整数。 定解区域内部的节点称为内点,记内点集R 为h 。边界与网格线的交点称为边界点,边界点全体记为h

33、 。与节点(,k j x y 沿x 方向或y 方向只差一个步长的点1(,k j x y 和1(,k j x y 称为节点(,k j x y 的相邻节点。如果一个内点的四个相邻节点均属于,称为正则内点,正内点的全体记为(1,至少有一个相邻节点不属于的内点称为非正则内点,非正则内点的全体记为(2。问题是要求出第一边值问题在全体内点上的数值解。 为简便,记:(,(,k j k j x y =,(,(,k j u k j u x y =,(,k j k j f f x y =。对正则内点(1(,k j ,由二阶中心差商公式:4422(4122(,22(4222(,(1,2(,(1,(,12(,12(,

34、(,1(,12k j x k j k j y k j u u k j u k j u k j h u x h y x h u u k j u k j u k j u x y y +-+-=-+-+-=-+ (16.5.2Poisson 方程2222(,u uf x y x y +=在点(,k j 处可表示为:,22(1,2(,(1,(,12(,(,1(,k j u k j u k j u k j u k j u k j u k j f R k j h +-+-+-+-+=+ (16.5.3 其中:4422(4(4122212(,(,(,1212(0,1k j k j x x h R k j u

35、 x h y u x y O h =+=+ (16.5.4 为其截断误差表示式,略去,(j k R ,即得与方程相近似的差分方程:1,1,1,1,2222k j k j k j kj k j kj k j u u u u u u f h +-+-+-+= (16.5.5式中方程的个数等于正则内点的个数,而未知数,k j u 则除了包含正则内点处解u 的近似值外,还包含一些非正则内点处u 的近似值,因而方程个数少于未知数个数。在非正则内点处Poisson 方程的差分近似不能按上式给出,需要利用边界条件得到。边界条件的处理可以有各种方案,下面介绍较简单的两种。(1 直接转移。用最接近非正则内点的边界点上的u 值作为该点上u 值的近似,这就是边界条件的直接转移。例如,点(,P k j 为非正则内点,其最接近的边界点为Q 点,则有: (2,(,k j u u Q Q k j = (16.3.6上式可以看作是用零次插值得到非正则内点处u 的近似值,容易求出,其截断误差为(+h O 。将上式代入,方程个数即与未知数个数相等。(2 线性插值。这种方案是通过用同一条网格线上与点P 相邻的边界点与内点作线性插值得到非正则内点(,P k j 处u 值的近似。由