《偏微分方程概述及運用matlab求解偏微分方程常見問題》

1、北京航空航天大學偏微分方程概述及運用matlab求解微分方程求解常見問題姓名徐敏學號57000211班級380911班2011年6月偏微分方程概述及運用matlab求解偏微分方程常見問題徐敏摘要偏微分方程簡介，matlab偏微分方程工具箱應用簡介，用這個工具箱解方程的過程是：確定待解的偏微分方程；確定邊界條件；確定方程所在域的幾何形狀；劃分有限元；解方程關鍵詞MATLAB偏微分方程程序如果一個微分方程中出現(xiàn)的未知函數(shù)只含有一個自變量，這個方程叫做常微分方程，也簡稱微分方程：如果一個微分方程中出現(xiàn)多元函數(shù)的偏導數(shù)，或者說如果未知函數(shù)和幾個變量有關，而且方程中出現(xiàn)未知函數(shù)對幾個變量的導數(shù)，那么這種

2、微分方程就是偏微分方程。一，偏微分方程概述偏微分方程是反映有關的未知變量關于時間的導數(shù)和關于空間變量的導數(shù)之間制約關系的等式。許多領域中的數(shù)學模型都可以用偏微分方程來描述，很多重要的物理、力學等學科的基本方程本身就是偏微分方程。早在微積分理論剛形成后不久，人們就開始用偏微分方程來描述、解釋或預見各種自然現(xiàn)象，并將所得到的研究方法和研究成果運用于各門科學和工程技術中，不斷地取得了顯著的成效，顯示了偏微分方程對于人類認識自然界基本規(guī)律的重要性。逐漸地，以物理、力學等各門科學中的實際問題為背景的偏微分方程的研究成為傳統(tǒng)應用數(shù)學中的一個最主要的內(nèi)容，它直接聯(lián)系著眾多自然現(xiàn)象和實際問題，不斷地提出和產(chǎn)生

3、出需要解決的新課題和新方法，不斷地促進著許多相關數(shù)學分支（如泛函分析、微分幾何、計算數(shù)學等）的發(fā)展，并從它們之中引進許多有力的解決問題的工具。偏微分方程已經(jīng)成為當代數(shù)學中的一個重要的組成部分，是純粹數(shù)學的許多分支和自然科學及工程技術等領域之間的一座重要的橋梁。在國外，對偏微分方程的應用發(fā)展是相當重視的。很多大學和研究單位都有應用偏微分方程的研究集體，并得到國家工業(yè)、科學部門及軍方、航空航天等方面的大力資助。比如在國際上有重大影響的美國的Courant研究所、法國的信息與自動化國立研究所等都集中了相當多的偏微分方程的研究人員，并把數(shù)學模型、數(shù)學方法、應用軟件及實際應用融為一體，在解決實際課題、推

4、動學科發(fā)展及加速培養(yǎng)人才等方面都起了很大的作用。在我國，偏微分方程的研究起步較晚。但解放后，在黨和國家的大力號召和積極支持下，我國偏微分方程的研究工作發(fā)展比較迅速，涌現(xiàn)出一批在這一領域中做出杰出工作的數(shù)學家，如谷超豪院士、李大潛院士等，并在一些研究方向上達到了國際先進水平。但總體來說，偏微分方程的研究隊伍的組織和水平、研究工作的廣度和深度與世界先進水平相比還有很大的差距。因此，我們必須繼續(xù)努力，大力加強應用偏微分方程的研究，逐步縮小與世界先進水平的差距二，偏微分方程的內(nèi)容偏微分方程是什么樣的？它包括哪些內(nèi)容？這里我們可從一個例子的研究加以介紹。弦振動是一種機械運動，當然機械運動的基本定律是質點

5、力學的F二ma，但是弦并不是質點，所以質點力學的定律并不適用在弦振動的研究上。然而，如果我們把弦細細地分成若干個極小極小的小段，每一小段抽象地看作是一個質點，這樣我們就可以應用質點力學的基本定律了。弦是指又細又長的彈性物質，比如弦樂器所用的弦就是細長的、柔軟的、帶有彈性的。演奏的時候，弦總是繃緊著具有一種張力，這種張力大于弦的重量幾萬倍。當演奏的人用薄片撥動或者用弓在弦上拉動，雖然只因其所接觸的一段弦振動，但是由于張力的作用，傳播到使整個弦振動起來。用微分的方法分析可得到弦上一點的位移是這一點所在的位置和時間為自變量的偏微分方程。偏方程又很多種類型，一般包括橢圓型偏微分方程、拋物型偏微分方程、

6、雙曲型偏微分方程。上述的例子是弦振動方程，它屬于數(shù)學物理方程中的波動方程，也就是雙曲型偏微分方程。偏微分方程的解一般有無窮多個，但是解決具體的物理問題的時候，必須從中選取所需要的解，因此，還必須知道附加條件。因為偏微分方程是同一類現(xiàn)象的共同規(guī)律的表示式，僅僅知道這種共同規(guī)律還不足以掌握和了解具體問題的特殊性，所以就物理現(xiàn)象來說，各個具體問題的特殊性就在于研究對象所處的特定條件，就是初始條件和邊界條件。拿上面所舉的弦振動的例子來說，對于同樣的弦的弦樂器，如果一種是以薄片撥動弦，另一種是以弓在弦上拉動，那么它們發(fā)出的聲音是不同的。原因就是由于“撥動”或“拉動”的那個“初始”時刻的振動情況不同，因此

7、產(chǎn)生后來的振動情況也就不同。天文學中也有類似情況，如果要通過計算預言天體的運動，必須要知道這些天體的質量，同時除了牛頓定律的一般公式外，還必須知道我們所研究的天體系統(tǒng)的初始狀態(tài)，就是在某個起始時間，這些天體的分布以及它們的速度。在解決任何數(shù)學物理方程的時候，總會有類似的附加條件。就弦振動來說，弦振動方程只表示弦的內(nèi)點的力學規(guī)律，對弦的端點就不成立，所以在弦的兩端必須給出邊界條件，也就是考慮研究對象所處的邊界上的物理狀況。邊界條件也叫做邊值問題。當然，客觀實際中也還是有“沒有初始條件的問題”，如定場問題（靜電場、穩(wěn)定濃度分布、穩(wěn)定溫度分布等），也有“沒有邊界條件的問題”，如著重研究不靠近兩端的那

8、段弦，就抽象的成為無邊界的弦了。在數(shù)學上，初始條件和邊界條件叫做定解條件。偏微分方程本身是表達同一類物理現(xiàn)象的共性，是作為解決問題的依據(jù);定解條件卻反映出具體問題的個性，它提出了問題的具體情況。方程和定解條件合而為一體，就叫做定解問題。求偏微分方程的定解問題可以先求出它的通解，然后再用定解條件確定出函數(shù)。但是一般來說，在實際中通解是不容易求出的，用定解條件確定函數(shù)更是比較困難的。偏微分方程的解法還可以用分離系數(shù)法，也叫做傅立葉級數(shù)；還可以用分離變數(shù)法，也叫做傅立葉變換或傅立葉積分。分離系數(shù)法可以求解有界空間中的定解問題，分離變數(shù)法可以求解無界空間的定解問題；也可以用拉普拉斯變換法去求解一維空間

9、的數(shù)學物理方程的定解。對方程實行拉普拉斯變換可以轉化成常微分方程，而且初始條件也一并考慮到，解出常微分方程后進行反演就可以了。應該指出，偏微分方程的定解雖然有以上各種解法，但是我們不能忽視由于某些原因有許多定解問題是不能嚴格解出的，只可以用近似方法求出滿足實際需要的近似程度的近似解。常用的方法有變分法和有限差分法。變分法是把定解問題轉化成變分問題，再求變分問題的近似解；有限差分法是把定解問題轉化成代數(shù)方程，然后用計算機進行計算；還有一種更有意義的模擬法，它用另一個物理的問題實驗研究來代替所研究某個物理問題的定解。雖然物理現(xiàn)象本質不同，但是抽象地表示在數(shù)學上是同一個定解問題，如研究某個不規(guī)則形狀

10、的物體里的穩(wěn)定溫度分布問題，在數(shù)學上是拉普拉斯方程的邊值問題，由于求解比較困難，可作相應的靜電場或穩(wěn)恒電流場實驗研究，測定場中各處的電勢，從而也解決了所研究的穩(wěn)定溫度場中的溫度分布問題。隨著物理科學所研究的現(xiàn)象在廣度和深度兩方面的擴展，偏微分方程的應用范圍更廣泛。從數(shù)學自身的角度看，偏微分方程的求解促使數(shù)學在函數(shù)論、變分法、級數(shù)展開、常微分方程、代數(shù)、微分幾何等各方面進行發(fā)展。從這個角度說，偏微分方程變成了數(shù)學的中心。三，用matlab解偏微分方程解偏微分方程不是一件輕松的事，但是偏微分方程在自然科學和工程領域中應用很廣，因此，我們可以運用matlab這個軟件來解決一些常見的偏微分方程。（一）

11、Matlab偏微分方程工具箱簡介。1，概述。本文只給出該工具箱的函數(shù)列表2，偏微分方程算法函數(shù)列表。adaptmesh生成自適應網(wǎng)絡及偏微分方程的解assemb生成邊界質量和剛度矩陣assema生成積分區(qū)域上質量和剛度矩陣assempde組成偏微分方程的剛度矩陣及右邊hyperbolic求解雙曲線型偏微分方程parabolic求解拋物線型偏微分方程pdeeig求解特征型偏微分方程pdenonlin求解非線性型微分方程poisolv利用矩陣格式快速求解泊松方程3，圖形界面函數(shù)。pdecirc畫圓pdeellip畫橢圓pdemdlcv轉化為版本1.0式的*.m文件pdepoly畫多邊形pderec

12、t畫矩形pdetool偏微分方程工具箱的圖形用戶界面4，幾何處理函數(shù)。csgchk檢查幾何矩陣的有效性csgdel刪除接近邊界的小區(qū)decsg將固定的幾何區(qū)域分解為最小區(qū)域initmesh產(chǎn)生最初的三角形網(wǎng)絡jigglemesh微調區(qū)域內(nèi)的三角形網(wǎng)絡poimesh在矩形區(qū)域上產(chǎn)生規(guī)則的網(wǎng)絡refinemesh細化三角形網(wǎng)絡wbound寫一個邊界描述文件wgeom寫一個幾何描述文件pdecont畫輪廓圖pdemesh畫偏微分方程的三角形網(wǎng)絡pdeplot畫偏微分方程的三角形網(wǎng)絡其中，A又稱為Laplace算子。這樣橢圓型偏微分方程可以簡單地寫pdesurf畫表面圖命令5，通用函數(shù)。pdetriq

13、三角形單元的品性度量poiasma邊界點對快速求解泊松方程的“貢獻”矩陣poicalc規(guī)范化的矩陣格式的點索引(二)Matlab偏微分方程工具箱應用。可以用詞工具箱求解如橢圓方程，雙曲線方程，特征值方程，拋物線方程。橢圓型偏微分方程的一般形式為一div(cVu)+au=f(x,t)其中：若u=u(x，x,，x,t)=u(x，r),Vu為u的梯度，則其定義為12nVu=散度div(v)的定義為66ddx，6x，6x1266+、6x6x12這樣，div(cVu)可以更明確地表示為div(v)=6)+6x丿ndiv(cVu)=6uc66x6x丿6x1126uc6x,26+U6x6x62+6x2丿nu

14、=cAu若c為常數(shù)，則進一步化簡為TOC o 1-5 h z6262div(cVu)=cI+I6x26x212特征值型偏微分方程c廠d2d2d2+dx2丿n+dx2dx212拋物型偏微分方程u+au=f(x,t)拋物型偏微分方程的一般形式為dud一div(cVu)+au=f(x,t)dt根據(jù)上面敘述，若c為常數(shù)，則該方程可以更簡單地寫為7dud一cdtd2d2d2+dx2dx2dx2)12nu+au=f(x,t)雙曲型偏微分方程雙曲型偏微分方程的一般形式為d2ud一div(cVu)+au=f(x,t)dt2若c為常數(shù)，則可以將該方程簡化為Td2u(d2d2d2d一cI+dt2Idx2dx2dx

15、2/12n三類方程的直接的區(qū)別在于u對t的導數(shù)的階次。u+au=f(x,t)若對t沒有求導，可以理解為其值為常數(shù)，故稱為橢圓型的。若取u對時間t的一階導數(shù)，則與u對x的二階導數(shù)直接構成了拋物線關系，故稱為拋物型偏微分方程。若取u對時間t的二階導數(shù)，稱其為雙曲型偏微分方程。特征值型偏微分方程為一div(cVu)+au-xdu對常數(shù)c該方程可以化簡為(d2d2d2)一cI+u+auxf(x,t)(dx2dx2dx2丿12n該方程是橢圓型偏微分方程的特例。pdetool的使用：在matlab命令窗口中鍵入pdetoo1窗口打開進入工作狀態(tài)，pdetool提供兩種解方程的方法，一種是通過函數(shù)，利用函數(shù)

16、可以編程也可以用命令行的方式解方程，兩一種是對pdetool窗口進行交互操作。一般來說，用函數(shù)解方程比較繁瑣，但是比較靈活：通過窗口交互操作比較簡單。解方程的全部過程以及結果都可以輸出保存為文本文件。限于文本的篇幅，我們主要介紹交互操作偏微分方程的方法。1.確定待解的偏微分方程。利用函數(shù)assempde可以對待解的偏微分方程加以描述。在交互操作中，為了方便用戶，pdetool把常見問題歸結為及各類型，可以再pdetool窗口的工具欄上找到選擇類型的彈出菜單，這些類型如下：通用問題通用系統(tǒng)(二維的偏微分方程組)結構力學：平面應力結構力學：平面應變靜電學靜磁學交流電電磁學直流電導電介質熱傳導擴散確

17、定問題類型后，可以再PDESpecification對話框填入c,a,f,d等系數(shù)，這樣就確定了待解的偏微分方程。2.確定邊界條件用函數(shù)assemb可以描述邊界條件。用pdetool提供的邊界條件對話框，在對話框里填入g,h,q,r等邊界條件。3.確定偏微分方程所在Q域的幾何圖形平面上波得散射問題。按照上面所說的解方程的過程，首先確定帶解的偏微分方程。散射是介質對入射波的反射。假定介質是均勻的，那么入射波在介質中傳播的速度是一個常數(shù)c?？梢杂煤瘮?shù)畫出0域的幾何圖形。Pdecirc:畫圓：pdeellip:畫橢圓：pderect:畫矩形：pdepoly：畫多邊形。無論哪種畫法，圖形一經(jīng)畫出，pd

18、etool就為這個圖形自動取名，并把代表圖形的名字放入Setformula窗口，在這個窗口，可以通過+，-圖形的名字現(xiàn)在對圖形的拓撲運算，以便構造復雜的0域幾何圖形4.劃分有限元對q域進行有限元的劃分函數(shù)有，initmesh:基本劃分；jigglemesh:采用jiggle方法劃分;refinemesh:精細劃分。在pdetool窗口中直接點擊劃分有限元的按鈕劃分有限元，劃分的方法與上面的函數(shù)想對應。5.解方程經(jīng)過1.11.4步驟后就可以解方程。解方程的函數(shù)有，adaptunesh:解方程的通用函數(shù)；poicalc：矩形有限元快速解橢圓形方程；poisolv:矩形有限元解橢圓形方程；parab

19、olic:解拋物線型方程：hyperbolic:解雙曲線型方程。在pdetool窗口中直接點擊解方程的按鈕即可解方程，解方程所耗費的時間在于有限元劃分的多少。(三)實例求解偏微分方程的邊值問題一、MATLAB支持的偏微分方程類型考慮平面有界區(qū)域D上的二階橢圓型PDE邊值問題：_VcVu)+au=f(0.1)其中(a力、(1)，一(2)a,f是D上的已知函數(shù)c是標量或2x2的函數(shù)方陣oxBy丿未知函數(shù)為U(x,y)(x,y)eD。它的邊界條件分為三類：1)Direchlet條件：hu=f2)Neumann條件：JcVu)+qu=g(3)混合邊界條件：在邊界qd上部分為Direchlet條件,為N

20、eumann條件。其中h,r,q,g,c是定義在邊界QD的已知函數(shù)，另外c也可以殳的函數(shù)矩陣，n是沿邊界的外法線的單位向量。在使用pdetool時要向它提供這些已知參數(shù)。二、例題例題1用pdetool求解-Au=1D:x2+y21u=0QD(0.2)(0.3)另外部分一個2*2(0.4)解：首先在MATLAB的工作命令行中鍵入pdetool，按回牟鍵確定于是出現(xiàn)PDEToolbox窗口，選GenenicScalar模式.(l)畫區(qū)域圓單擊橢圓工具按鈕，大致在(0，0)位置單擊鼠標右鍵，拖拉鼠標到適當位置松開。為了保證所繪制的圓是標準的單位園，在所繪園上雙擊，打開ObjectDialog對話框，

21、精確地輸入圓心坐標X-center為0、Y-center為0及半徑Radius為l，然后單擊OK按鈕，這樣單位畫已畫好(2)設置邊界條件單擊工具邊界模式按鈕，圖形邊界變紅，逐段雙擊邊界，打開Boundarycondition對話框輸入邊界條件對于同一類型的邊界，可以按Shift鍵，將多個邊界同時選擇，統(tǒng)一設邊界條件.本題選擇Dirichlet條件，輸入h為1,r為0。，然后單擊OK按鈕也可以單擊Boundary菜單中SpocifyBoundaryCondition選項，打開BoundaryCondition對話框輸入邊界條件(3)設置方程單擊偏微分方程按鈕，打開PDESpecification

22、對話框，選擇方程類型本題選Ellintic(橢圓型)，輸入c為1,a為0,f為1,然后單擊OK按鈕(4)網(wǎng)格剖分單擊網(wǎng)格工具，或者單擊Mesh菜單中InitializeMesh項，可進行初始網(wǎng)格剖分這時在PDEToolbox窗口下方的狀態(tài)欄內(nèi)顯示出初始網(wǎng)格的節(jié)點數(shù)和三角形單元數(shù)本題節(jié)點數(shù)為144個，三角形單元數(shù)為254個(圖?)。如果要細化網(wǎng)格，單擊細化工具，或者單擊Mesh菜單中RefineMesh選項，節(jié)點數(shù)成為541個，三角形單元數(shù)為1016個。(5)解方程單擊解方程工具，或者單擊Solve菜單中SolvePDE選項，可求得方程數(shù)值解并用彩色圖形顯示。單擊作圖工具，或者單擊Plot菜單中

23、Parameter選項，出現(xiàn)Plotselection對話框.從中選擇于Height(3-Dplot)，然后單擊Plot按鈕，方程的圖形解如圖?？所示。除了作定解問題解u的圖形外，也可以作|gradu|等圖形Info:Selectanewplot,orchangemodetoalterPDEZmesh,orboundaries.PDEToolbox-UntitledFileEditOptionsDrawBoundaryPDEMeshSolvePlotWindowHelp|田|O|方|OG|pDe|A|血J=|氣矗enericSc砌三11%：心Setformula:-1.5-1-0.50.51.

24、6）輸出網(wǎng)格節(jié)點的編號、單元編號以及節(jié)點坐標單擊Mesh菜單中ShowNodeLabels選項，再單擊網(wǎng)格工具，即可顯示節(jié)點編號（圖？）。若要輸出節(jié)點坐標，只需單擊Mesh菜單中ExportMesh選項，這時打開的Export對話框中的默認值為pet,這里p、e、t分別表示point（點）、edges（邊）、triangles（三角形）數(shù)據(jù)變量，單擊OK按鈕，然后在MATLAB命令行鍵入p,即可以顯示按節(jié)點編號排列的坐標;鍵入e再回車則顯示邊界數(shù)據(jù)矩陣（7維數(shù)組）；鍵入t按回車則顯示三角形單元數(shù)據(jù)矩陣（4維數(shù)組）。點、邊、單元的部分輸出為：p=Columns1through11-0.7071-0.70710.7071-0.7071-1.00000.00001.00000.00000.7071-0.7071-0.9808-0.9239-0.8315-0.0000-1.000001.00000.70710.7071-0.1951-0.3827-0.5556Columns1through111.00009.000010.000011.00005.000012.000013.000014.00002.00001