Ch7-求解常微分方程Word版

1、225-232第七章 求解偏微分方程常微分方程的求解是在高等數(shù)學(xué)里面講過的，所以大家比較熟悉，對于偏微分方程的求解可能就有點(diǎn)陌生了。事實(shí)上偏微分方程在工程上有著更廣泛地應(yīng)用，例如描述液體在多孔介質(zhì)中的擴(kuò)散，聲學(xué)和電磁學(xué)中諧波的傳播，熱在固體中的傳導(dǎo)等過程的微分方程都是偏微分方程。所以求解偏微分方程在工程實(shí)際上有著非常重要的價(jià)值。手工求解常微分方程就已經(jīng)非常麻煩和復(fù)雜了，那么求解偏微分方程就更加的復(fù)雜和繁冗了。所幸的是有MATLAB這個(gè)計(jì)算工具來幫我們解決這一麻煩問題，MATLAB中有一個(gè)專門用來求解片微分方程的工具箱PDE工具箱。這個(gè)工具箱不但提供有豐富的命令函數(shù)，使求解偏微分方程便的簡單靈活

2、，而且該提供了一個(gè)求解偏微分方程的圖形用戶界面系統(tǒng)（GUI），使得整個(gè)求解過程更加的人性化。特別是對于初學(xué)者來說GUIjiang更容易被接受，操作也更方便，所以偏微分方程的求解這一部分我們將主要介紹PDE工具箱的GUI系統(tǒng)。7-1 偏微分方程的特點(diǎn)對于n階常微分方程我們知道它的解取決于n個(gè)任意常數(shù)，例如一階常微分方程的解可以表示為：,c為任意常數(shù)。但是對于偏微分方程來說就不是由多少個(gè)常數(shù)來確定的了，例如偏微分方程的解可表示為： 其中w(x)和v(y)為兩個(gè)連續(xù)的任意可微函數(shù)。我們看到偏微分方程的解可能是非常多的，與常微分方程的解依賴于若干任意常數(shù)相比，它的自由度要大的多，對于多維偏微分方程的解

3、更是這樣。正是由于這個(gè)原因，一般來說偏微分方程的解很難用通式表達(dá)出來。事實(shí)上我們常用到的是偏微分方程在某種特定條件下的解，這樣靠著這些特定條件的約束我們就可以把偏微分方程的解表示出來了。我們把這些幫助確定偏微分方程特解的條件叫做定解條件，由于自變量在多維空間中的變化，其變化的區(qū)域非常復(fù)雜，所以在區(qū)域邊界上給出的定解條件也更加形式多樣，我們一般稱給定在區(qū)域邊界上的定解條件為邊界條件。我們看到邊界條件歲于求解偏微分方程非常有用。由于偏微分方程的求解一般來說太過復(fù)雜，所以現(xiàn)在沒有一個(gè)對所有的偏微分進(jìn)行求解的理論，所謂的求解偏微分方程也只是對于某些人們比較熟悉的類型的偏微分方程進(jìn)行求解。不同類型的偏微

4、分方程代表不同的物理現(xiàn)象，例如方程表示內(nèi)部沒有熱源的情況下物體的溫度分布應(yīng)滿足的條件，t表示時(shí)間變量，x、y、z分別表示物體的三維空間坐標(biāo)，這里1 / 11，K為物體的熱傳導(dǎo)系數(shù)，Q為該物體的熱的容量。再比如方程表示聲波或者振動(dòng)在介質(zhì)中的傳播，同樣t表示振動(dòng)傳播的時(shí)間變量，x、y、z表示介質(zhì)的三維空間坐標(biāo)，而a則表示振動(dòng)在該介質(zhì)中傳播的速度。人們通過對偏微分方程的研究，把常見的偏微分方程歸結(jié)為四種類型，即：橢圓型方程、拋物型方程和混合型方程。前面講過的熱傳導(dǎo)方程是一個(gè)拋物型方程，波動(dòng)方程則是一個(gè)雙曲型方程。如果熱傳導(dǎo)方程不考慮時(shí)間因素，則方程變?yōu)椋@個(gè)方程就是一個(gè)橢圓方程了，混合型方程顧名思義

5、就是一個(gè)方程具有上面三種類型方程中至少兩種的特點(diǎn)的方程了。上面說過，邊界條件對于求解偏微分方程很關(guān)鍵，在求解過程中最常用的地定解條件有兩種，即：狄利克雷條件（又叫第一類邊界條件）、諾曼條件（又叫第二類邊界條件）。狄利克雷條件常寫為，即滿足給定區(qū)域邊界的待求函數(shù)u為已知函數(shù)g；對于熱傳導(dǎo)方程狄利克雷條件的具體意義就是，已知物體的表面溫度分布函數(shù)為g ;對于波動(dòng)方程則狄利克雷條件的具體意義就是，波在介質(zhì)的表面（或介質(zhì)的給定區(qū)域的表面）的振動(dòng)函數(shù)u為已知函數(shù)g。諾曼條件常寫為，為已知函數(shù)，是關(guān)于的外法線導(dǎo)數(shù)。對于熱傳導(dǎo)方程諾曼條件的具體意義為，沿給定區(qū)域邊界的外法向的溫度變化率為已知函數(shù)；對于波動(dòng)方

6、程則諾曼條件的具體意義為，沿給定區(qū)域邊界的外法向的振動(dòng)速率為已知函數(shù)。7-2 PDE toolbox求解問題的背景知識(shí)1求解方程的類型使用PDE工具箱求解偏微分方程的基本類型有橢圓型方程、雙曲型方程、拋物型方程、特征值方程、橢圓型方程組和非線形橢圓型方程組。PDE工具箱對這些方程的公式表述為：橢圓型方程：雙曲型方程：拋物型方程：特征值方程：這里是微分算子，表示對多元函數(shù)求全微分運(yùn)算。是給定的平面有界區(qū)域，參數(shù)c、a、f和待求未知數(shù)u都是定義在上的函數(shù)。參數(shù)d是定義在上的復(fù)函數(shù)，是待求的特征值。在拋物型方程和雙曲型方程中，參數(shù)函數(shù)c、a、f、d可以是時(shí)間t的函數(shù)?？梢允褂梅蔷€形解題器求解非線形橢

7、圓型方程：這里c、a和f都是未知的待求函數(shù)u的函數(shù)。另外，PDE工具箱提供的解題器都可以處理下面的方程組：利用命令還可以求解高階方程組。2邊界條件我們知道常用的兩個(gè)邊界條件狄利克雷條件和諾曼條件，PDE工具箱對它們的公式表示為：Dirichlet條件：Neumann條件為：這里是上的單位外法向矢量，h、r、q和g定義在上的復(fù)值函數(shù)（對于特征值問題g=0,r=0）。對于非線性問題，系數(shù)g、q、h和 r可以依賴于u；對于拋物型方程和雙曲型方程系數(shù)還可以依賴于時(shí)間t。對于方程組情形，狄利克雷邊界條件表示為：一般諾曼邊界條件表示為：混合型邊界條件表示為：這里計(jì)算時(shí)要注意滿足狄利克雷條件。7-3 GUI

8、求解偏微分方程的過程7-3-1 GUI求解問題的一般程序在命令窗口里輸入pdetool，回車，PDE工具箱的圖形用戶界面（GUI）系統(tǒng)就啟動(dòng)了。從定義一個(gè)偏微分方程問題到完成解偏微分方程的定解，整個(gè)過程大致可以分為六個(gè)階段，它們是：（1） 繪制平面有界區(qū)域，被MATLAB稱之為“Draw模式”。系統(tǒng)提供了一系列的實(shí)體模型，包括矩形、圓、橢圓和多邊形，我們可以通過公式把這些實(shí)體模型組合起來，生成我們需要的平面區(qū)域。（2） 定義邊界，被MATLAB稱之為“Boundary模式”。聲明不同邊界段的邊界條件。（3） 定義偏微分方程，被MATLAB稱之為“PDE模式”。確定方程的類型和方程的系數(shù)c、a、

9、f和d。根據(jù)具體情況（例如材料的性質(zhì)）,我們還可以在不同的子區(qū)域上獨(dú)立的聲明不同的系數(shù)。（4） 網(wǎng)格化區(qū)域，被MATLAB稱之為“Mesh模式”。我們可以控制自動(dòng)生成網(wǎng)格的參數(shù)，還可以對生成的網(wǎng)格進(jìn)行對次的細(xì)化，使網(wǎng)格分割的更細(xì)更合理。（5） 解偏微分方程，被MATLAB稱之為“Solve模式”。對于橢圓型方程，我們可以激活并控制非線性自適應(yīng)解題器來處理非線性方程；對于拋物線方程和雙曲線方程，設(shè)置初始邊界條件后可以求出給定時(shí)刻t的解；對于特征值問題，可以求出給定區(qū)間上的特征值。求解完成后我們還可以返回到第四步，對網(wǎng)絡(luò)進(jìn)一步細(xì)化，進(jìn)行再次求解。（6）計(jì)算結(jié)果的可視化。我我們可以通過設(shè)置系統(tǒng)提供的

10、對話框，顯示所求的解的表面圖、網(wǎng)絡(luò)圖、等高線圖和箭頭梯度圖（quiver）。對于拋物型和雙曲型問題還可以進(jìn)行動(dòng)畫演示。如何通過上面的六個(gè)階段來具體的求解一個(gè)偏微分方程，我們將在下一節(jié)中通過一個(gè)示例給出具體的操作步驟。7-3-2 求解問題的示例讓我們從求解最簡單的偏微分方程泊松方程開始，泊松方程的表達(dá)式為：這里是拉普拉斯算子，注意和微分算子區(qū)分開。定解條件我們采用狄利克雷邊界條件。我們看到泊松方程屬于橢圓型偏微分方程，它可以表示內(nèi)部沒有熱源的物體的熱的分布，還可以表示靜止的電磁場的分布。首先在MATLAB的命令窗口輸入pdetool，回車確認(rèn)，這樣圖形用戶界面系統(tǒng)（GUI）就啟動(dòng)了。接下來我們將

11、按照上一節(jié)所說的六個(gè)步驟開始操作了，啟動(dòng)后的GUI窗口如圖7-1所示。圖7-1 PDE工具箱的GUI 窗口在打開的窗口中，點(diǎn)擊菜單Options，選中Grid選項(xiàng)，于是窗口中便出現(xiàn)了網(wǎng)格，這樣方便我們在繪圖時(shí)確定幾何圖形的位置和大小。用鼠標(biāo)左鍵單擊工具欄里的畫圓按鈕（左邊第四個(gè)），使其處于選中狀態(tài)（表現(xiàn)為下沉），然后在窗口的繪圖區(qū)域按下左鍵拖動(dòng)鼠標(biāo)，于是一個(gè)橢圓就出現(xiàn)了，系統(tǒng)給它命名為E1。再選中工具欄里的畫矩形的按鈕（左邊第二個(gè)），在繪圖區(qū)域的中心位置拖動(dòng)鼠標(biāo)，于是一個(gè)矩形就出現(xiàn)了，系統(tǒng)給它命名為R1。把光標(biāo)移到E1上面，雙擊鼠標(biāo)左鍵，這時(shí)會(huì)彈出一個(gè)對話框，用來定義E1的位置和大小，設(shè)置如圖

12、7-2所示。圖7-2 幾何實(shí)體E1屬性對話框然后選擇OK按鈕確認(rèn)。用同樣的辦法調(diào)出R1的屬性對話框，設(shè)置如下，然后確認(rèn)，如圖7-3所示。圖7-3 R1屬性對話框接下來把窗口中Set formula欄中的公式修改為E1-R1，這樣幾何實(shí)體的繪制就完成了，畫完以后窗口如圖7-4所示。我們用鼠標(biāo)點(diǎn)住畫好的幾何實(shí)體，然后拖動(dòng)鼠標(biāo)就可以移動(dòng)實(shí)體的位置，另外還可以把畫好的幾何實(shí)體用m文件的形式保存起來。圖7-4 繪好的幾何實(shí)體接下來設(shè)置邊界。用鼠標(biāo)左鍵單擊工具欄里的按鈕，這樣邊界就生成了，如圖7-5所示。圖7-5 生成邊界條件我們看到邊界是分段的并且有方向（帶有箭頭）。把光標(biāo)移到一段邊界上，然后雙擊鼠標(biāo)左

13、鍵，就會(huì)出現(xiàn)一個(gè)對話框，讓我們設(shè)置該段的邊界條件，如圖7-6所示。本例中所有的邊界都采用Dirichlet條件的默認(rèn)值。圖7-6 邊界條件對話框我們既可以采用上面的辦法一段一段的設(shè)置邊界條件，也可以一次設(shè)置好邊界條件相同的邊界段。先按下shift鍵，然后用鼠標(biāo)左鍵分別單擊條件性同的邊段，這樣這些邊界段都被選中了，再使用上面的辦法彈出對話框設(shè)置參數(shù)，于是被選中的段都被設(shè)置了。選擇方程的類型。單擊工具欄里的PDE 按鈕，打開方程 類型對話框，由于 泊松方程適宜中橢圓型方程，所以在對話框中選擇Elliptic選項(xiàng)。方程的參數(shù)設(shè)置如圖7-7所示，單擊OK按鈕確認(rèn)。圖7-7 方程類型對話框?qū)υO(shè)定的區(qū)域進(jìn)行網(wǎng)格剖分。單擊工具欄里的三角號按鈕，就可以進(jìn)行網(wǎng)格剖分了。我們還可以單擊工具欄里網(wǎng)格剖分按鈕右邊的按鈕（一個(gè)三角號里面還有一個(gè)小三角的按鈕）對網(wǎng)格剖分進(jìn)一步細(xì)化，網(wǎng)格細(xì)化后的GUI窗口如圖7-8所示。圖7-8 網(wǎng)格細(xì)化后的圖像求解方程。單擊工具欄里的等號按鈕，開始求解。這時(shí)窗口里就顯示出了泊松方程在給定邊界條件下的圖形解，如圖7-9所示。圖7-9 泊松方程的圖形解該圖中的圖像是通過顏色的深淺來表示z坐標(biāo)的大小的。我們不但能用平面圖來表示方程的解，還可以使用多種方式可視化方程的解。例如在圖7