第二十章 偏微分方程的數(shù)值解

1、第二十章 偏微分方程的數(shù)值解自然科學與工程技術中種種運動發(fā)展過程與平衡現(xiàn)象各自遵守一定的規(guī)律。這些規(guī)律的定量表述一般地呈現(xiàn)為關于含有未知函數(shù)及其導數(shù)的方程。我們將只含有未知多元函數(shù)及其偏導數(shù)的方程，稱之為偏微分方程。方程中出現(xiàn)的未知函數(shù)偏導數(shù)的最高階數(shù)稱為偏微分方程的階。如果方程中對于未知函數(shù)和它的所有偏導數(shù)都是線性的，這樣的方程稱為線性偏微分方程，否則稱它為非線性偏微分方程。初始條件和邊界條件稱為定解條件，未附加定解條件的偏微分方程稱為泛定方程。對于一個具體的問題，定解條件與泛定方程總是同時提出。定解條件與泛定方程作為一個整體，稱為定解問題。§1 偏微分方程的定解問題各種物理性質(zhì)的

2、定常（即不隨時間變化）過程，都可用橢圓型方程來描述。其最典型、最簡單的形式是泊松(Poisson)方程 （1） 特別地，當時，即為拉普拉斯(Laplace)方程，又稱為調(diào)和方程 （2）帶有穩(wěn)定熱源或內(nèi)部無熱源的穩(wěn)定溫度場的溫度分布，不可壓縮流體的穩(wěn)定無旋流動及靜電場的電勢等均滿足這類方程。Poisson方程的第一邊值問題為 （3）其中為以為邊界的有界區(qū)域，為分段光滑曲線，稱為定解區(qū)域，分別為上的已知連續(xù)函數(shù)。第二類和第三類邊界條件可統(tǒng)一表示成 （4）其中為邊界的外法線方向。當時為第二類邊界條件，時為第三類邊界條件。在研究熱傳導過程，氣體擴散現(xiàn)象及電磁場的傳播等隨時間變化的非定常物理問題時，常常

3、會遇到拋物型方程。其最簡單的形式為一維熱傳導方程 （5）方程（5）可以有兩種不同類型的定解問題：初值問題（也稱為Cauchy問題） （6）初邊值問題 （7）其中為已知函數(shù)，且滿足連接條件 問題（7）中的邊界條件稱為第一類邊界條件。第二類和第三類邊界條件為 （8）其中。當時，為第二類邊界條件，否則稱為第三類邊界條件。雙曲型方程的最簡單形式為一階雙曲型方程 （9）物理中常見的一維振動與波動問題可用二階波動方程 （10）描述，它是雙曲型方程的典型形式。方程（10）的初值問題為 （11）邊界條件一般也有三類，最簡單的初邊值問題為如果偏微分方程定解問題的解存在，唯一且連續(xù)依賴于定解數(shù)據(jù)（即出現(xiàn)在方程和定

4、解條件中的已知函數(shù)），則此定解問題是適定的。可以證明，上面所舉各種定解問題都是適定的。§2 偏微分方程的差分解法差分方法又稱為有限差分方法或網(wǎng)格法，是求偏微分方程定解問題的數(shù)值解中應用最廣泛的方法之一。它的基本思想是：先對求解區(qū)域作網(wǎng)格剖分，將自變量的連續(xù)變化區(qū)域用有限離散點（網(wǎng)格點）集代替；將問題中出現(xiàn)的連續(xù)變量的函數(shù)用定義在網(wǎng)格點上離散變量的函數(shù)代替；通過用網(wǎng)格點上函數(shù)的差商代替導數(shù)，將含連續(xù)變量的偏微分方程定解問題化成只含有限個未知數(shù)的代數(shù)方程組（稱為差分格式）。如果差分格式有解，且當網(wǎng)格無限變小時其解收斂于原微分方程定解問題的解，則差分格式的解就作為原問題的近似解（數(shù)值解）。

5、因此，用差分方法求偏微分方程定解問題一般需要解決以下問題：（i）選取網(wǎng)格；（ii）對微分方程及定解條件選擇差分近似，列出差分格式；（iii）求解差分格式；（iv）討論差分格式解對于微分方程解的收斂性及誤差估計。下面我們只對偏微分方程的差分解法作一簡要的介紹。2.1 橢圓型方程第一邊值問題的差分解法以Poisson方程（1）為基本模型討論第一邊值問題的差分方法。考慮Poisson方程的第一邊值問題（3） 取分別為方向和方向的步長，以兩族平行線 將定解區(qū)域剖分成矩形網(wǎng)格。節(jié)點的全體記為。定解區(qū)域內(nèi)部的節(jié)點稱為內(nèi)點，記內(nèi)點集為。邊界與網(wǎng)格線的交點稱為邊界點，邊界點全體記為。與節(jié)點沿方向或方向只差一個

6、步長的點和稱為節(jié)點的相鄰節(jié)點。如果一個內(nèi)點的四個相鄰節(jié)點均屬于，稱為正則內(nèi)點，正則內(nèi)點的全體記為，至少有一個相鄰節(jié)點不屬于的內(nèi)點稱為非正則內(nèi)點，非正則內(nèi)點的全體記為。我們的問題是要求出問題（3）在全體內(nèi)點上的數(shù)值解。為簡便記，記。對正則內(nèi)點，由二階中心差商公式 Poisson方程（1）在點處可表示為 （12）在式（12）中略去，即得與方程（1）相近似的差分方程 (13) 式（13）中方程的個數(shù)等于正則內(nèi)點的個數(shù)，而未知數(shù)則除了包含正則內(nèi)點處解的近似值，還包含一些非正則內(nèi)點處的近似值，因而方程個數(shù)少于未知數(shù)個數(shù)。在非正則內(nèi)點處Poisson方程的差分近似不能按式（13）給出，需要利用邊界條件得到

7、。邊界條件的處理可以有各種方案，下面介紹較簡單的兩種。(i) 直接轉移(ii) 線性插值由式（13）所給出的差分格式稱為五點菱形格式，實際計算時經(jīng)常取，此時五點菱形格式可化為 (14)簡記為 （15）其中。 求解差分方程組最常用的方法是同步迭代法，同步迭代法是最簡單的迭代方式。除邊界節(jié)點外，區(qū)域內(nèi)節(jié)點的初始值是任意取定的。例1 用五點菱形格式求解Laplace方程第一邊值問題 其中。取。當時，利用點構造的差分格式 （16）稱為五點矩形格式，簡記為 （17） 其中。2.2 拋物型方程的差分解法以一維熱傳導方程（5） 為基本模型討論適用于拋物型方程定解問題的幾種差分格式。首先對平面進行網(wǎng)格剖分。分

8、別取為方向與方向的步長，用兩族平行直線，將平面剖分成矩形網(wǎng)格，節(jié)點為。為簡便起見，記，。2.2.1 微分方程的差分近似在網(wǎng)格內(nèi)點處，對分別采用向前、向后及中心差商公式，對采用二階中心差商公式，一維熱傳導方程（5）可分別表示為由此得到一維熱傳導方程的不同的差分近似 （18） （19） （20）初、邊值條件的處理為用差分方程求解定解問題（6），（7）等，還需對定解條件進行離散化。對初始條件及第一類邊界條件，可直接得到 （21） （22）其中。對第二、三類邊界條件則需用差商近似。下面介紹兩種較簡單的處理方法。（i）在左邊界處用向前差商近似偏導數(shù)，在右邊界處用向后差商近似偏導數(shù)，即 即得邊界條件（8）

9、的差分近似為 （23）（ii）用中心差商近似，即 則得邊界條件的差分近似為 （24）這樣處理邊界條件，誤差的階數(shù)提高了，但式（24）中出現(xiàn)定解區(qū)域外的節(jié)點和，這就需要將解拓展到定解區(qū)域外?？梢酝ㄟ^用內(nèi)節(jié)點上的值插值求出和，也可以假定熱傳導方程（5）在邊界上也成立，將差分方程擴展到邊界節(jié)點上，由此消去和。幾種常用的差分格式下面我們以熱傳導方程的初邊值問題（7）為例給出幾種常用的差分格式。(i) 古典顯式格式為便于計算，令，式（18）改寫成以下形式將式（18）與（21），（22）結合，我們得到求解問題（7）的一種差分格式 （25）由于第0層上節(jié)點處的值已知，由式（25）即可算出在第一層上節(jié)點處的近

10、似值。重復使用式（25），可以逐層計算出各層節(jié)點的近似值。（ii）古典隱式格式將（19）整理并與式（21），（22）聯(lián)立，得差分格式如下（26）其中。雖然第0層上的值仍為已知，但不能由式（30）直接計算以上各層節(jié)點上的值，故差分格式（26）稱為古典隱式格式。（iii）杜福特弗蘭克爾（DoFortFrankel）格式DoFortFrankel格式是三層顯式格式，它是由式（24）與（25），（26）結合得到的。具體形式如下：(27)用這種格式求解時，除了第0層上的值由初值條件（21）得到，必須先用二層格式求出第1層上的值，然后再按格式（27）逐層計算。2.3 雙曲型方程的差分解法對二階波動方程（1

11、0） 如果令，則方程（10）可化成一階線性雙曲型方程組 （28）記，則方程組（28）可表成矩陣形式 （29）矩陣有兩個不同的特征值，故存在非奇異矩陣，使得 作變換，方程組（29）可化成 （30）方程組（30）由兩個獨立的一階雙曲型方程聯(lián)立而成。因此下面主要討論一階雙曲型方程的差分解法?？紤]一階雙曲型方程的初值問題 (31)與拋物型方程的討論類似，仍將平面剖分成矩形網(wǎng)格。取方向步長為，方向步長為，網(wǎng)格線為。為簡便起見，記，。 以不同的差商近似偏導數(shù)，可以得到方程（9）的不同的差分近似 （32） （33） （34）結合離散化的初始條件，可以得到幾種簡單的差分格式。§3 Matlab解法M

12、atlab中的偏微分方程（PDE）工具箱是用有限元法尋求典型偏微分方程的數(shù)值近似解，該工具箱求解偏微分方程具體步驟與用有限元方法求解偏微分方程的過程是一致的，包括幾個步驟，即幾何描述、邊界條件描述、偏微分方程類型選擇、有限元劃分計算網(wǎng)格、初始化條件輸入，最后給出偏微分方程的數(shù)值解（包括畫圖）。下面我們討論的方程是定義在平面上的有界區(qū)域上，區(qū)域的邊界記作。3.1 方程類型Matlab工具箱可以解決下列類型的偏微分方程：（i）橢圓型偏微分方程 其中和未知的可以是上的復值函數(shù)。（ii）拋物型偏微分方程 其中可以依賴于時間。（iii）雙曲型偏微分方程 （iv）特征值問題 其中是未知的特征值，是上的復值

13、函數(shù)。（v）非線性橢圓偏微分方程 其中可以是的函數(shù)。（vi）方程組 3.2 邊界條件邊界條件有如下三種：（i）Dirichlet條件： on 。（ii）Neumann條件： on 。這里為區(qū)域的單位外法線，是定義在上的復值函數(shù)。對于二維方程組情形，Dirichlet邊界條件為 ， ；Neumann邊界條件為： （iii）對于偏微分方程組，混合邊界條件為 這里的計算是使得滿足Dirichlet邊界條件。3.3 求解偏微分方程例2 求解泊松方程 ，求解區(qū)域為單位圓盤，邊界條件為在圓盤邊界上。解 它的精確解為 下面求它的數(shù)值解，編寫程序如下：%(1)問題定義g='circleg' %

14、單位圓b='circleb1' %邊界上為零條件c=1;a=0;f=1;%（2）產(chǎn)生初始的三角形網(wǎng)格p,e,t=initmesh(g); %（3）迭代直至得到誤差允許范圍內(nèi)的合格解error=; err=1;while err > 0.01,p,e,t=refinemesh(g,p,e,t);u=assempde(b,p,e,t,c,a,f); %求得數(shù)值解exact=(1-p(1,:).2-p(2,:).2)/4;err=norm(u-exact',inf);error=error err;end%結果顯示subplot(2,2,1),pdemesh(p,e,t

15、);subplot(2,2,2),pdesurf(p,t,u)subplot(2,2,3),pdesurf(p,t,u-exact')例3 考慮最小表面問題 在 ，在圓盤邊界上。解 這是橢圓型方程，其中，編寫程序如下：g='circleg'b='circleb2'c='1./sqrt(1+ux.2+uy.2)'rtol=1e-3;p,e,t=initmesh(g);p,e,t=refinemesh(g,p,e,t);u=pdenonlin(b,p,e,t,c,0,0,'Tol',rtol);pdesurf(p,t,u)例4

16、 求解正方形區(qū)域上的熱傳導方程 初始條件為邊界條件為Dirichlet條件。 解 這里是拋物型方程，其中。編寫程序如下：%(1)問題定義g='squareg' %定義正方形區(qū)域b='squareb1' %邊界上為零條件c=1;a=0;f=0;d=1;%（2）產(chǎn)生初始的三角形網(wǎng)格p,e,t=initmesh(g); %(3)定義初始條件u0=zeros(size(p,2),1);ix=find(sqrt(p(1,:).2+p(2,:).2)<0.4);u0(ix)=1%(4)在時間段為0到0.1的20個點上求解nframe=20;tlist=linspace

17、(0,0.1,nframe);u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);%(5)動畫圖示結果for j=1:nframe pdesurf(p,t,u1(:,j); mv(j)=getframe;endmovie(mv,10) 例5 求解正方形區(qū)域上的波方程 初始條件為，邊界條件為在上滿足Dirichlet條件，在上滿足Neumann條件。 解 這里是雙曲型方程，其中。編寫程序如下：%(1)問題定義g='squareg' %定義正方形區(qū)域b='squareb3' %定義邊界c=1;a=0;f=0;d=1;%（2）產(chǎn)生初始的三角形

18、網(wǎng)格p,e,t=initmesh(g); %(3)定義初始條件x=p(1,:)'y=p(2,:)'u0=atan(cos(pi*x);ut0=3*sin(pi*x).*exp(cos(pi*y);%(4)在時間段為0到5的31個點上求解n=31;tlist=linspace(0,5,n);uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);%(5)動畫圖示結果for j=1:n pdesurf(p,t,uu(:,j); mv(j)=getframe;endmovie(mv,10) 例6 求解泊松方程 ，求解區(qū)域為單位圓盤，邊界條件為在圓盤邊界

19、上。解 它的精確解為 。下面求它的數(shù)值解，編寫程序如下：g='circleg' b='circleb1' c=1;a=0;f='circlef'p,e,t=initmesh(g);p,e,t=refinemesh(g,p,e,t);u=assempde(b,p,e,t,c,a,f); exact=-1/(2*pi)*log(sqrt(p(1,:).2+p(2,:).2);subplot(2,2,1),pdemesh(p,e,t);subplot(2,2,2),pdesurf(p,t,u)subplot(2,2,3),pdesurf(p,t,u-e

20、xact')3.4 偏微分方程的pdetool解法對于一般的區(qū)域，任意邊界條件的偏微分方程，我們可以利用Matlab中pdetool提供的偏微分方程用戶圖形界面解法。pdetool提供的用戶圖形界面解法的使用步驟如下：（i）在Matlab命令窗口運行pdetool，出現(xiàn)PDE Toolbox界面。（ii）用鼠標點一下工具欄上的“PDE"按鈕，在彈出的對話框中定義偏微分方程。（iii）用鼠標點一下工具欄上的區(qū)域按鈕，在下面的坐標系中畫出偏微分方程的大致定解區(qū)域。（iv）雙擊（iii）中畫出的大致區(qū)域，在彈出的對話框中精確定位定解區(qū)域。（v）用鼠標點一下工具欄上的邊界按鈕“”，畫出區(qū)域的邊界。（vi）雙擊坐標系中的區(qū)域邊界，定義偏微分方程的邊界條件。（vii）用鼠標點工具欄上的剖分按鈕，對求解區(qū)域進行剖分。（viii）如果求拋物型或雙曲型方程的數(shù)值解，還需要通過“solve”菜單下的“parameters”選項設置初值條件。（ix）用鼠標點一下工具欄上的“”按鈕，就畫出偏微分方程數(shù)值解的圖形。通過“solve”菜單下的“Export Solution”選項可以把數(shù)值解u輸出到Matlab的工作間。（x）如要畫出數(shù)值解的三維圖形，需要設置“plot