Ch7-求解常微分方程Word版

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

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

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

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

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

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

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

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

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

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

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

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

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