偏微分方程的數(shù)值解法

1、第十六章 偏微分方程的數(shù)值解法科學研究和工程技術中的許多問題可建立偏微分方程的數(shù)學模型。包含多個自變量的微分方程稱為偏微分方程(partial differential equation,簡稱PDE 。偏微分方程問題,其求解是十分困難的。除少數(shù)特殊情況外,絕大多數(shù)情況均難以求出精確解。因此,近似解法就顯得更為重要。本章僅介紹求解各類典型偏微分方程定解問題的差分方法。16.1 幾類偏微分方程的定解問題一個偏微分方程的表示通常如下:(,x x x y y y x y A B C f x y += (16.1.1 式中,A B C 是常數(shù),稱為擬線性(quasilinear數(shù)。通常,存在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選取網(wǎng)格: 圖16.2.1 差分示意圖首先對定解區(qū)域(,0D x t x t =-+作網(wǎng)格剖分,最簡單常用一種網(wǎng)格是用兩族分別平行于x 軸與t 軸的等距直線:k x x k h =,(0,1,2,0,1,2,j t t j k j = (16.2.2將D 分成許多小矩形區(qū)域。這些直線稱為網(wǎng)格線,其交點稱為網(wǎng)格點,也稱為節(jié)點,h 和分別稱作x 方向和t 方向的步長。這種網(wǎng)格稱為矩形網(wǎng)格。(1 對微分方程及定解條件選擇差分近似,列出差分格式。如果用向前差商表示一階偏

4、導數(shù),即: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 ,用差分格式求其數(shù)值解,(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 比較知,按前一個差分格式所求得的數(shù)值解不收斂到初值問題的解,而后一個差分格式的解收斂到原問題的解。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 方程組由二個獨立的一階雙曲型方程聯(lián)立而成。因此本節(jié)主要討論一階雙曲型方程的差分解法。下面給出如下波動方程和邊界條件的差

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 -平面進行網(wǎng)格剖分。分別取,h 為x 方向與t 方向的步長,用兩族平行直線(0,1,2,k x x kh k =,(0,1,2j t t j j =,將x t -平面剖分成矩形網(wǎng)格,節(jié)點為:(,(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 微分方程的差分近似在網(wǎng)格內(nèi)點(,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上述差分方程所用到的節(jié)點各不相同。其截斷誤差分別為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差分方程用到四個節(jié)點。由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 =時,才成立。 綜上可知,用不同的差商公式可以得到微分方程的不同的差分近似。構造差分格式的關鍵在于使其具有相容性、收斂性和穩(wěn)定性。前面三個方程都具有相容性,而此方程則要在一定條件下才具有相容性。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 =處用向前差商近似偏導數(shù)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 。誤差的階數(shù)提高了,但出現(xiàn)定解區(qū)域外的節(jié)點(1,j -和(1,n j +,這就需要將解拓展到定解區(qū)域外?？梢酝ㄟ^用

15、內(nèi)節(jié)點上的u 值插值求出1,j u -和1,n j u +,也可以假定熱傳導方程在邊界上也成立,將差分方程擴展到邊界節(jié)點上,由此消去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、=上節(jié)點處的u 值已知(0,k k u =,由此即可算出u 在第一層(1j =上節(jié)點處的近似值,1k u 。重復使用此式,可以逐層計算出所有的,k j u ,因此此差分格式稱為典顯式格式。又因式中只出現(xiàn)相鄰兩個時間層的節(jié)點,故此式是二層顯式格式。(2 古典隱式格式 將式:,11,1,220k j k j k j k j k ju u u u u a h-+-+-= (16.4.21 整理并與初始條件及第一類邊界條件式聯(lián)立,得差分格式如下:,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 值仍為已知,但不能由上式直接計算以上各層節(jié)點上的值,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、式。此方程組是三對角方程組,且系數(shù)矩陣嚴格對角占優(yōu),故解存在唯一。(3 Richardson 格式 Richardson 格式是將式:,1,11,1,220k j k jk jk j k ju u uu u ah +-+-+-=整理后與初始條件及第一類邊界條件式聯(lián)立。其計算公式為:,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這種差分格式中所涉及的節(jié)點出現(xiàn)在1,1+-j j j 三層

20、上,故為三層顯式格式。Richardson 格式是一種完全不穩(wěn)定的差分格式,因此它在實際計算中是不能采用的。(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 +,將上式與初始條件及第一類邊界條件式聯(lián)立并整理,得差分格式:,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 此格式涉及到六個節(jié)點,它又是隱式格式,故稱為六點隱式格式。與古典隱式格式類似,用六點格式由第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 =-。此方程組的系數(shù)矩陣嚴格對角占優(yōu),故仍可用追趕法求解。例16.4.1 用古典顯式格式求初邊值問題。222003,03(,003(0,0,(3,903u u t x t x u x x x u t u t t -=的數(shù)值解,取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 =如此逐層計算,得全部節(jié)點上的數(shù)值解為: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 ,求解區(qū)間(,:0,0R x t x a t b =內(nèi)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 ,求解區(qū)間(,:0,0R x t x a t b =內(nèi)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、分解法本節(jié)以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 =將定解區(qū)域剖分成矩形網(wǎng)格。節(jié)點的全體記為: (,k j k j R x y x k h y j k j =為整數(shù)。 定解區(qū)域內(nèi)部的節(jié)點稱為內(nèi)點,記內(nèi)點集R 為h 。邊界與網(wǎng)格線的交點稱為邊界點,邊界點全體記為h

33、 。與節(jié)點(,k j x y 沿x 方向或y 方向只差一個步長的點1(,k j x y 和1(,k j x y 稱為節(jié)點(,k j x y 的相鄰節(jié)點。如果一個內(nèi)點的四個相鄰節(jié)點均屬于,稱為正則內(nèi)點,正內(nèi)點的全體記為(1,至少有一個相鄰節(jié)點不屬于的內(nèi)點稱為非正則內(nèi)點,非正則內(nèi)點的全體記為(2。問題是要求出第一邊值問題在全體內(nèi)點上的數(shù)值解。 為簡便,記:(,(,k j k j x y =,(,(,k j u k j u x y =,(,k j k j f f x y =。對正則內(nèi)點(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式中方程的個數(shù)等于正則內(nèi)點的個數(shù),而未知數(shù),k j u 則除了包含正則內(nèi)點處解u 的近似值外,還包含一些非正則內(nèi)點處u 的近似值,因而方程個數(shù)少于未知數(shù)個數(shù)。在非正則內(nèi)點處Poisson 方程的差分近似不能按上式給出,需要利用邊界條件得到。邊界條件的處理可以有各種方案,下面介紹較簡單的兩種。(1 直接轉移。用最接近非正則內(nèi)點的邊界點上的u 值作為該點上u 值的近似,這就是邊界條件的直接轉移。例如,點(,P k j 為非正則內(nèi)點,其最接近的邊界點為Q 點,則有: (2,(,k j u u Q Q k j = (16.3.6上式可以看作是用零次插值得到非正則內(nèi)點處u 的近似值,容易求出,其截斷誤差為(+h O 。將上式代入,方程個數(shù)即與未知數(shù)個數(shù)相等。(2 線性插值。這種方案是通過用同一條網(wǎng)格線上與點P 相鄰的邊界點與內(nèi)點作線性插值得到非正則內(nèi)點(,P k j 處u 值的近似。由