Matlab解微分方程

1、第十六章 偏微分方程的數(shù)值解法科學(xué)研究和工程技術(shù)中的許多問題可建立偏微分方程的數(shù)學(xué)模型。包含多個自變量的微分方程稱為偏微分方程(partial differential equation)，簡稱PDE。偏微分方程問題，其求解是十分困難的。除少數(shù)特殊情況外，絕大多數(shù)情況均難以求出精確解。因此，近似解法就顯得更為重要。本章僅介紹求解各類典型偏微分方程定解問題的差分方法。16.1 幾類偏微分方程的定解問題一個偏微分方程的表示通常如下： (16.1.1)式中，是常數(shù)，稱為擬線性(quasilinear)數(shù)。通常，存在3種擬線性方程： 雙曲型(hyperbolic)方程：； 拋物線型(parabolic

2、)方程：； 橢圓型(ellliptic)方程：。16.1.2 雙曲型方程最簡單形式為一階雙曲型方程： (16.1.2)物理中常見的一維振動與波動問題可用二階波動方程： (16.1.3)描述，它是雙曲型方程的典型形式。方程的初值問題為： (16.1.4)邊界條件一般有三類，最簡單的初邊值問題為： (16.1.5)16.1.3 拋物型方程其最簡單的形式為一維熱傳導(dǎo)方程： (16.1.8)方程可以有兩種不同類型的定解問題：(1) 初值問題： (16.1.6)(2) 初邊值問題： (16.1.7)其中，為已知函數(shù)，且滿足連接條件： (16.1.8)邊界條件為第一類邊界條件。第二類和第三類邊界條件為：

3、(16.1.9)其中，。當(dāng)時，為第二類邊界條件，否則稱為第三類邊界條件。16.1.4 橢圓型方程其最典型、最簡單的形式是泊松(Poisson)方程 (16.1.10)特別地，當(dāng)時，即為拉普拉斯(Laplace)方程，又稱為調(diào)和方程： (16.1.11)Poisson方程的第一邊值問題為： (16.1.12)其中為以為邊界的有界區(qū)域，為分段光滑曲線，稱為定解區(qū)域，分別為，上的已知連續(xù)函數(shù)。第二類和第三類邊界條件可統(tǒng)一表示為： (16.1.13)其中為邊界的外法線方向。當(dāng)時為第二類邊界條件，時為第三類邊界條件。16.2 差分方法的基本概念差分方法又稱為有限差分方法或網(wǎng)格法，是求偏微分方程定解問題的

4、數(shù)值解中應(yīng)用最廣泛的方法之一。它的基本思想是：先對求解區(qū)域作網(wǎng)格剖分，將自變量的連續(xù)變化區(qū)域用有限離散點(diǎn)(網(wǎng)格點(diǎn))集代替；將問題中出現(xiàn)的連續(xù)變量的函數(shù)用定義在網(wǎng)格點(diǎn)上離散變量的函數(shù)代替；通過用網(wǎng)格點(diǎn)上函數(shù)的差商代替導(dǎo)數(shù)，將含連續(xù)變量的偏微分方程定解問題化成只含有限個未知數(shù)的代數(shù)方程組(稱為差分格式)。如果差分格式有解，且當(dāng)網(wǎng)格無限變小時其解收斂于原微分方程定解問題的解，則差分格式的解就作為原問題的近似解(數(shù)值解)。因此，用差分方法求偏微分方程定解問題一般需要解決以下問題：(1) 選取網(wǎng)格；(2) 對微分方程及定解條件選擇差分近似，列出差分格式；(3) 求解差分格式；(4) 討論差分格式解對于微

5、分方程解的收斂性及誤差估計(jì)。下面，用一個簡單的例子來說明用差分方法求解偏微分方程問題的一般過程及差分方法的基本概念。設(shè)有一階雙曲型方程初值問題。 (16.2.1)選取網(wǎng)格：圖16.2.1 差分示意圖首先對定解區(qū)域作網(wǎng)格剖分，最簡單常用一種網(wǎng)格是用兩族分別平行于軸與軸的等距直線： , (16.2.2)將分成許多小矩形區(qū)域。這些直線稱為網(wǎng)格線，其交點(diǎn)稱為網(wǎng)格點(diǎn)，也稱為節(jié)點(diǎn)，和分別稱作方向和方向的步長。這種網(wǎng)格稱為矩形網(wǎng)格。(1) 對微分方程及定解條件選擇差分近似，列出差分格式。如果用向前差商表示一階偏導(dǎo)數(shù)，即： (16.2.3)其中。方程： (16.2.4)在節(jié)點(diǎn)處可表示為： (16.2.5)其中

6、：。由于當(dāng)足夠小時，在式中略去，就得到一個與方程相近似的差分方程： (16.2.6)此處，可看作是問題的解在節(jié)點(diǎn)處的近似值。同初值條件： (16.2.7)結(jié)合，就得到求問題的數(shù)值解的差分格式。式： (16.2.8)稱為差分方程的截?cái)嗾`差。如果一個差分方程的截?cái)嗾`差為，則稱差分方程對是階精度，對是階精度的。顯然，截?cái)嗾`差的階數(shù)越大，差分方程對微分方程的逼近越好。若網(wǎng)格步長趨于0時，差分方程的截?cái)嗾`差也趨于0，則稱差分方程與相應(yīng)的微分方程是相容的。這是用差分方法求解偏微分方程問題的必要條件。如果當(dāng)網(wǎng)格步長趨于0時，差分格式的解收斂到相應(yīng)微分方程定解問題的解，則稱這種差分格式是收斂的。16.3 雙曲

7、型方程的差分解法16.3.1 一階雙曲型方程的差分格式考慮一階雙曲型方程的初值問題： (16.3.1)將平面剖分成矩形網(wǎng)格，取方向步長為方向步長為，網(wǎng)格線為： , 為簡便，記： , 以不同的差商近似偏導(dǎo)數(shù)，可以得到方程的不同的差分近似 (16.3.2) (16.3.3) (16.3.4)截?cái)嗾`差分別為，與。結(jié)合離散化的初始條件，可以得到幾種簡單的差分格式： (16.3.5) (16.3.6) (16.3.7)其中：。如果已知第層節(jié)點(diǎn)上的值，按上面三種格式就可求出第層上的值。因此，這三種格式都是顯式格式。如果對采用向后差商，采用向前差商，則方程可化成： (16.3.8)相應(yīng)的差分格式為： (16

8、.3.9)此差分格式是一種隱式格式，必須通過解方程組才能由第層節(jié)點(diǎn)上的值求出第層節(jié)點(diǎn)上的值。例16.3.1 對初值問題： 其中：，用差分格式求其數(shù)值解，取。解：記，由初始條件： 按差分格式： 計(jì)算公式為：。計(jì)算結(jié)果略。如果用差分格式： 求解，計(jì)算公式為： 計(jì)算結(jié)果略。與準(zhǔn)確解： 比較知，按前一個差分格式所求得的數(shù)值解不收斂到初值問題的解，而后一個差分格式的解收斂到原問題的解。16.3.2 波動方程的差分格式對二階波動方程： (16.3.1)如果令，則方程可化成一階線性雙曲型方程組： (16.3.2)記，則方程組可表成矩陣形式： (16.3.3)矩陣有兩個不同的特征值，故存在非奇異矩陣，使得：

9、(16.3.4)作變換，方程組可化為： (16.3.5)方程組由二個獨(dú)立的一階雙曲型方程聯(lián)立而成。因此本節(jié)主要討論一階雙曲型方程的差分解法。 下面給出如下波動方程和邊界條件的差分格式： (16.3.6) (16.3.7) 將矩形劃分成個小矩形，長寬分別為：，形成一個網(wǎng)格，如圖16.3.1所示。圖16.3.1 網(wǎng)格圖可通過差分方程的方法計(jì)算近似值： 在連續(xù)行內(nèi)，網(wǎng)格點(diǎn)的真實(shí)值為。 求和的中心差分為： (16.3.8) (16.3.9)在每一行的網(wǎng)格間距是均勻的：，而且。同時，它在每一列也是均勻的，。接下來，將和去掉，用(16.3.8)和(16.3.9)中的近似，并按順序代入(16.3.6)式，可

10、得到差分方程： (16.3.10)可用它來近似(16.3.6)式。為方便起見，可將代入上式，可得： (16.3.11)設(shè)行和的近似值已知，可用上式求網(wǎng)格的行： (16.3.12)式中，。根據(jù)上式左邊的4個已知值可得到近似值，如圖16.3.2所示。圖16.3.2 波動方程的空格樣板 用上式時，必須注意，如果計(jì)算的某個階段帶來的誤差最終會越來越小，則方法是穩(wěn)定的。為了保證上式的穩(wěn)定性，必然使。還存在其他一些差分方程方法以，稱為隱格式法，它們更難實(shí)現(xiàn)，但對無限制。16.3.3 差分方法求解波動方程的MATLAB程序 求解區(qū)間，以(16.3.7)為邊界條件的波動方程的差分方法程序。*function

11、U = finedif(f,g,a,b,c,n,m) %Input - f=u(x,0) as a string 'f'% - g=ut(x,0) as a string 'g'% - a and b right endpoints of 0,a and 0,b% - c the constant in the wave equation% - n and m number of grid points over 0,a and 0,b%Output - U solution matrix; analogous to Table 10.1%Initialize

12、parameters and U h = a/(n-1);k = b/(m-1);r = c*k/h;r2=r2;r22=r2/2;s1 = 1 - r2;s2 = 2 - 2*r2;U = zeros(n,m); %Comput first and second rows for i=2:n-1 U(i,1)=feval(f,h*(i-1); U(i,2)=s1*feval(f,h*(i-1)+k*feval(g,h*(i-1) . +r22*(feval(f,h*i)+feval(f,h*(i-2);end %Compute remaining rows of U for j=3:m, f

13、or i=2:(n-1), U(i,j) = s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1)-U(i,j-2); endend U=U'*16.4 拋物型方程的差分解法以一維熱傳導(dǎo)方程： (16.4.1)為基本模型討論適用于拋物型方程定解問題的幾種差分格式。16.4.1 差分格式的建立首先對平面進(jìn)行網(wǎng)格剖分。分別取為方向與方向的步長，用兩族平行直線，將平面剖分成矩形網(wǎng)格，節(jié)點(diǎn)為：。為簡便，記： ，。16.4.2 微分方程的差分近似在網(wǎng)格內(nèi)點(diǎn)處，對分別采用向前、向后及中心差商公式： (16.4.2)一維熱傳導(dǎo)方程可分別表示為： (16.4.3) (16.4.4

14、) (16.4.5)由此得到一維熱傳導(dǎo)方程的不同差分近似： (16.4.6) (16.4.7) (16.4.8)上述差分方程所用到的節(jié)點(diǎn)各不相同。其截?cái)嗾`差分別為，和。因此，它們都與一維熱傳導(dǎo)方程相容。如果將式：中的用代替，則可得到又一種差分近似： (16.4.9)差分方程用到四個節(jié)點(diǎn)。由Taylor公式容易得出： (16.4.10)故其的截?cái)嗾`差為。因而不是對任意的，此差分方程都能逼近熱傳導(dǎo)方程：，僅當(dāng)時，才成立。綜上可知，用不同的差商公式可以得到微分方程的不同的差分近似。構(gòu)造差分格式的關(guān)鍵在于使其具有相容性、收斂性和穩(wěn)定性。前面三個方程都具有相容性，而此方程則要在一定條件下才具有相容性。1

15、6.4.3 初、邊值條件的處理為用差分方法求解定解問題初值問題： (16.4.11)初邊值問題： (16.4.12)還需對定解條件進(jìn)行離散化。對初始條件及第一類邊界條件，可直接得到： ， ， 式中：。對第二、三類邊界條件： (16.4.13)需用差分近似。下面介紹兩種較簡單的處理方法。在左邊界處用向前差商近似偏導(dǎo)數(shù)，在右邊界處用向后差商近似，即： ， (16.4.14)則得邊界條件的差分近似為： ， (16.4.15)其截?cái)嗾`差為。(2) 用中心差商近似,即： ， (16.4.16)則得邊界條件的差分近似為： ， (16.4.17)其截?cái)嗾`差為。誤差的階數(shù)提高了，但出現(xiàn)定解區(qū)域外的節(jié)點(diǎn)和，這就

16、需要將解拓展到定解區(qū)域外?？梢酝ㄟ^用內(nèi)節(jié)點(diǎn)上的值插值求出和，也可以假定熱傳導(dǎo)方程在邊界上也成立，將差分方程擴(kuò)展到邊界節(jié)點(diǎn)上，由此消去和。16.4.4 幾種常用的差分格式以熱傳導(dǎo)方程的初邊值問題： (16.4.18)為例給出幾種常用的差分格式。(1) 古典顯式格式令，則： (16.4.19)可改寫成： 。將其與初始條件及第一類邊界條件： ， ，結(jié)合，我們得到求解此問題的一種差分格式： (16.4.20)由于第0層上節(jié)點(diǎn)處的值已知，由此即可算出在第一層上節(jié)點(diǎn)處的近似值。重復(fù)使用此式，可以逐層計(jì)算出所有的，因此此差分格式稱為典顯式格式。又因式中只出現(xiàn)相鄰兩個時間層的節(jié)點(diǎn)，故此式是二層顯式格式。(2)

17、 古典隱式格式將式： (16.4.21)整理并與初始條件及第一類邊界條件式聯(lián)立，得差分格式如下： (16.4.22)其中：。雖然第0層上的值仍為已知，但不能由上式直接計(jì)算以上各層節(jié)點(diǎn)上的值，必須通過解下列線性方程組： (16.4.23)式中，。 才能由計(jì)算，故此差分格式稱為古典隱式格式。此方程組是三對角方程組，且系數(shù)矩陣嚴(yán)格對角占優(yōu)，故解存在唯一。(3) Richardson格式Richardson格式是將式： 整理后與初始條件及第一類邊界條件式聯(lián)立。其計(jì)算公式為： 6.4.24)這種差分格式中所涉及的節(jié)點(diǎn)出現(xiàn)在三層上，故為三層顯式格式。Richardson格式是一種完全不穩(wěn)定的差分格式，因此

18、它在實(shí)際計(jì)算中是不能采用的。(4) 杜福特-弗蘭克爾 (DoFort-Frankel) 格式DoFort-Frankel格式也是三層顯式格式，它是由式： 與初始條件及第一類邊界條件式結(jié)合得到的。具體形式如下： (16.4.25)用這種格式求解時，除了第0層上的值由初值條件得到，必須先用二層格式求出第1層上的值，然后再按上式逐層計(jì)算。(5) 六點(diǎn)隱式格式對二階中心差商公式： 如果用在與處的二階中心差商的平均值近似在處的值，即： 同時在點(diǎn)處的值也用中心差商： 近似，即： 這樣又得到熱傳導(dǎo)方程的一種差分近似： (16.4.26)其截?cái)嗾`差為，將上式與初始條件及第一類邊界條件式聯(lián)立并整理，得差分格式：

19、 (14.5.27)此格式涉及到六個節(jié)點(diǎn)，它又是隱式格式，故稱為六點(diǎn)隱式格式。與古典隱式格式類似，用六點(diǎn)格式由第層的值計(jì)算第層的值時，需求解三對角方程組： (16.4.28)式中，。此方程組的系數(shù)矩陣嚴(yán)格對角占優(yōu)，故仍可用追趕法求解。例16.4.1 用古典顯式格式求初邊值問題。 的數(shù)值解，取。解：這里：，。由格式： 可得到： 將初值代入上式，即可算出：將邊界條件，及上述結(jié)果代入又可求得：如此逐層計(jì)算，得全部節(jié)點(diǎn)上的數(shù)值解為：16.4.5 前向差分法求解熱傳導(dǎo)方程的MATLAB程序 設(shè)，其中，而且，其中，求解區(qū)間內(nèi)的近似解。*function U=forwdif(f,c1,c2,a,b,c,n,

20、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.4%Initialize parameters and U h=a/(n-1);k=b/(m-

21、1);r=c2*k/h2;s=1-2*r;U=zeros(n,m); %Boundary conditions U(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 U for j=2:m for i=2:n-1 U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1); endend U=U'*16.4.6 Crank-Nicholson求解熱傳導(dǎo)方程的MATLAB程序 設(shè)，其中，而且，其中，求

22、解區(qū)間內(nèi)的近似解。*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% - n and m number of grid points over 0,a and 0,b%Output - U solution matrix; analogous to Table 10.5%I

23、nitialize parameters and U h=a/(n-1);k=b/(m-1);r=c2*k/h2;s1=2+2/r;s2=2/r-2;U=zeros(n,m); %Boundary conditions U(1,1:m)=c1;U(n,1:m)=c2; %Generate first row U(2:n-1,1)=feval(f,h:h:(n-2)*h)' %Form the diagonal and off-diagonal elemnts of A and %the constant vector B and solve tridiagonal system AX=

24、B Vd(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:m for i=2:n-1 Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1); end X=trisys(Va,Vd,Vc,Vb); U(1:n,j)=X'end U=U'*16.5 橢圓型方程第一邊值問題的差分解法本節(jié)以Poisson方程為基本模型討論第一邊值問題的差分方法。16.5.1 差分格式的建立考慮Poiss

25、on方程的第一邊值問題： (16.5.1)圖16.5.1取分別為方向和方向的步長，如圖所示，以兩族平行線：，將定解區(qū)域剖分成矩形網(wǎng)格。節(jié)點(diǎn)的全體記為： 。定解區(qū)域內(nèi)部的節(jié)點(diǎn)稱為內(nèi)點(diǎn)，記內(nèi)點(diǎn)集為。邊界與網(wǎng)格線的交點(diǎn)稱為邊界點(diǎn)，邊界點(diǎn)全體記為。與節(jié)點(diǎn)沿方向或方向只差一個步長的點(diǎn)和稱為節(jié)點(diǎn)的相鄰節(jié)點(diǎn)。如果一個內(nèi)點(diǎn)的四個相鄰節(jié)點(diǎn)均屬于，稱為正則內(nèi)點(diǎn)，正內(nèi)點(diǎn)的全體記為，至少有一個相鄰節(jié)點(diǎn)不屬于的內(nèi)點(diǎn)稱為非正則內(nèi)點(diǎn)，非正則內(nèi)點(diǎn)的全體記為。問題是要求出第一邊值問題在全體內(nèi)點(diǎn)上的數(shù)值解。 為簡便，記： ，。對正則內(nèi)點(diǎn)，由二階中心差商公式： (16.5.2)Poisson方程在點(diǎn)處可表示為： (16.5.3)其

26、中： (16.5.4) 為其截?cái)嗾`差表示式，略去，即得與方程相近似的差分方程： (16.5.5)式中方程的個數(shù)等于正則內(nèi)點(diǎn)的個數(shù)，而未知數(shù)則除了包含正則內(nèi)點(diǎn)處解的近似值外，還包含一些非正則內(nèi)點(diǎn)處的近似值，因而方程個數(shù)少于未知數(shù)個數(shù)。在非正則內(nèi)點(diǎn)處Poisson方程的差分近似不能按上式給出，需要利用邊界條件得到。邊界條件的處理可以有各種方案，下面介紹較簡單的兩種。(1) 直接轉(zhuǎn)移。用最接近非正則內(nèi)點(diǎn)的邊界點(diǎn)上的值作為該點(diǎn)上值的近似，這就是邊界條件的直接轉(zhuǎn)移。例如，點(diǎn)為非正則內(nèi)點(diǎn)，其最接近的邊界點(diǎn)為點(diǎn)，則有： (16.3.6)上式可以看作是用零次插值得到非正則內(nèi)點(diǎn)處的近似值，容易求出，其截?cái)嗾`差為

27、。將上式代入，方程個數(shù)即與未知數(shù)個數(shù)相等。(2) 線性插值。這種方案是通過用同一條網(wǎng)格線上與點(diǎn)相鄰的邊界點(diǎn)與內(nèi)點(diǎn)作線性插值得到非正則內(nèi)點(diǎn)處值的近似。由點(diǎn)與的線性插值確定的近似值，得： (16.5.7)其中，其截?cái)嗾`差為。將其與方程相近似的差分方程聯(lián)立，得到方程個數(shù)與未知數(shù)個數(shù)相等的方程組，求解此方程組可得Poisson方程第一邊值問題的數(shù)值解。上面所給出的差分格式稱為五點(diǎn)菱形格式， (16.5.8)實(shí)際計(jì)算時經(jīng)常取，此時五點(diǎn)菱形格式可化為： (16.5.9)簡記為： (16.5.10)其中：。例16.5.1 用五點(diǎn)菱形格式求解拉普拉斯(Laplace)方程第一邊值問題。 其中：。取。解：網(wǎng)格中

28、有四個內(nèi)點(diǎn)，均為正則內(nèi)點(diǎn)。由五點(diǎn)菱形格式，得方程組： 代入邊界條件： 其解為：，當(dāng)時，對利用點(diǎn)，構(gòu)造的差分格式，稱為五點(diǎn)矩形格式，簡記為： (16.5.11)其中：，其截?cái)嗾`差為： (16.5.12)五點(diǎn)菱形格式與矩形格式的截?cái)嗾`差均為，稱它們具有二階精度。如果用更多的點(diǎn)構(gòu)造差分格式，其截?cái)嗾`差的階數(shù)可以提高，如利用菱形格式及矩形格式所涉及的所有節(jié)點(diǎn)構(gòu)造出的九點(diǎn)格式就是具有四階精度的差分格式。16.5.2 用Dirichlet法求解Laplace方程的MATLAB程序 求解區(qū)間內(nèi)的近似解。而且滿足條件：，其中，而且，其中，而且存在整數(shù)和，使得：。*function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1) %Input - f1,f2,