偏微分方程的MATLAB解法

PAGEPAGE34引言偏微分方程定解問題有著廣泛的應(yīng)用背景。人們用偏微分方程來描述、解釋或者預見各種自然現(xiàn)象，并用于科學和工程技術(shù)的各個領(lǐng)域fll。然而，對于廣大應(yīng)用工作者來說，從偏微分方程模型出發(fā)，使用有限元法或有限差分法求解都要耗費很大的工作量，才能得到數(shù)值解?，F(xiàn)在，MATLABPDEToolbox已實現(xiàn)對于空間二維問題高速、準確的求解過程。偏微分方程如果一個微分方程中出現(xiàn)的未知函數(shù)只含一個自變量，這個方程叫做常微分方程，也簡稱微分方程；如果一個微分方程中出現(xiàn)多元函數(shù)的偏導數(shù)，或者說如果未知函數(shù)和幾個變量有關(guān)，而且方程中出現(xiàn)未知函數(shù)對幾個變量的導數(shù)，那么這種微分方程就是偏微分方程。常用的方法有變分法和有限差分法。變分法是把定解問題轉(zhuǎn)化成變分問題，再求變分問題的近似解；有限差分法是把定解問題轉(zhuǎn)化成代數(shù)方程，然后用計算機進行計算；還有一種更有意義的模擬法，它用另一個物理的問題實驗研究來代替所研究某個物理問題的定解。雖然物理現(xiàn)象本質(zhì)不同，但是抽象地表示在數(shù)學上是同一個定解問題，如研究某個不規(guī)則形狀的物體里的穩(wěn)定溫度分布問題，由于求解比較困難，可作相應(yīng)的靜電場或穩(wěn)恒電流場實驗研究，測定場中各處的電勢，從而也解決了所研究的穩(wěn)定溫度場中的溫度分布問題。隨著物理科學所研究的現(xiàn)象在廣度和深度兩方面的擴展，偏微分方程的應(yīng)用范圍更廣泛。從數(shù)學自身的角度看，偏微分方程的求解促使數(shù)學在函數(shù)論、變分法、級數(shù)展開、常微分方程、代數(shù)、微分幾何等各方面進行發(fā)展。從這個角度說，偏微分方程變成了數(shù)學的中心。一、MATLAB方法簡介及應(yīng)用1.1MATLAB簡介MATLAB是美國MathWorks公司出品的商業(yè)數(shù)學軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級技術(shù)計算語言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分。1.2Matlab主要功能數(shù)值分析數(shù)值和符號計算工程與科學繪圖控制系統(tǒng)的設(shè)計與仿真數(shù)字圖像處理數(shù)字信號處理通訊系統(tǒng)設(shè)計與仿真財務(wù)與金融工程1.3優(yōu)勢特點1)高效的數(shù)值計算及符號計算功能，能使用戶從繁雜的數(shù)學運算分析中解脫出來；2)具有完備的圖形處理功能，實現(xiàn)計算結(jié)果和編程的可視化；3)友好的用戶界面及接近數(shù)學表達式的自然化語言，使學者易于學習和掌握；4)功能豐富的應(yīng)用工具箱(如信號處理工具箱、通信工具箱等)，為用戶提供了大量方便實用的處理工具。1.4MATLAB產(chǎn)品族可以用來進行以下各種工作●數(shù)值分析●數(shù)值和符號計算●工程與科學繪圖●控制系統(tǒng)的設(shè)計與仿真●數(shù)字圖像處理技術(shù)●數(shù)字信號處理技術(shù)●通訊系統(tǒng)設(shè)計與仿真●財務(wù)與金融工程●管理與調(diào)度優(yōu)化計算（運籌學）MATLAB的應(yīng)用范圍非常廣，包括信號和圖像處理、通訊MATLAB在通訊系統(tǒng)設(shè)計與仿真的應(yīng)用、控制系統(tǒng)設(shè)計、測試和測量、財務(wù)建模和分析以及計算生物學等眾多應(yīng)用領(lǐng)域。附加的工具箱（單獨提供的專用MATLAB函數(shù)集）擴展了MATLAB環(huán)境，以解決這些應(yīng)用領(lǐng)域內(nèi)特定類型的問題。二、Laplacian算子簡介Laplacian算子：Poisson方程(ellipptic)：Laplacian算子的特征值問題：Heatequation(parabolic):Waveequation(hyperbolic):五點離散:Poisson方程離散：特征值問題：熱方程：波動方程：波動方程：橢圓方程：特征值方程：熱方程：波動方程：熱方程：波動方程：三、Matlab解偏微分方程解偏微分方程不是一件輕松的事情，但是偏微分方程在自然科學和工程領(lǐng)域應(yīng)用很廣，因此，研究解偏微分方程的方法、開發(fā)解偏微分方程的工具是數(shù)學和計算機領(lǐng)域中的一項重要工作。MATLAB提供了專門用于解二維偏微分方程的工具箱，使用這個工具箱，一方面解偏微分方程，另一方面，可以讓我們學習如何把求解數(shù)學問題的過程與方法工程化。應(yīng)當承認，我們國家在數(shù)學軟件的開發(fā)方面還比較落后，MATLAB是當今世界上最好的數(shù)學軟件之一，通過對這個軟件的認識，有助于研發(fā)我們自己的數(shù)學軟件。MATLAB的偏微分方程工具箱名字叫pdetool，它采用有限元法解偏微分方程。用這個工具箱可以解如下方程。橢圓方程拋物線方程雙曲線方程特征值方程所有的方程都在二維平面Ω域上。方程中，▽是Laplace算子，u是待解的未知函數(shù)，c,a,f是已知的實值標量函數(shù)，d是已知的復值函數(shù)，λ是未知的特征值。MATLAB提供了兩種方法Ⅲ解決PDE問題，一是pdepe函數(shù)，它可以求解一般的PDEs，具有較大的通用性，但只支持命令行形式調(diào)用。二是PDE工具箱，可以求解特殊PDE問題，但有較大的局限性，比如只能求解二階PDE問題，并且不能解決偏微分方程組。它提供了GUI界面，可以從繁雜的編程中解脫出來，同時還可以通過FileSaveAs直接生成M代碼。3.1函數(shù)解法pdepe函數(shù)介紹它的調(diào)用格式為sol=pdepe(m，@pdefun，@pdeic，@pdebc,x,t)輸入?yún)?shù):@pdefun：是PDE的問題描述函數(shù)@pdebc：是PDE的邊界條件描述函數(shù)@pdeic：是PDE的初值條件輸出參數(shù)：sol：是一個三維數(shù)組，sol(:,:,i)表示ui的解，換句話說uk對應(yīng)x(i)和t(j)時的解為sol(i,j，k)，通過sol，我們可以使用pdeval()直接計算某個點的函數(shù)值。實例講解例：試求解下面的偏微分其中， ，且滿足初始條件及邊界條件解：（1）對照給出的偏微分方程，根據(jù)標注形式，則原方程可以改寫為可見m=0,且％目標PDE函數(shù)function[c，f，s]--pdefun(x，t，u，du)c=[1；1]；f_[0．024*du(1)；0．17*du(2)]；temp=u(1)一u(2)；s=[一1；1]．*(exp(5．73*temp)一exp(-11．46*temp))；（2）邊界條件改寫為下邊界上邊界％邊界條件函數(shù)function[pa．qa，pb，qb]--pdebe(xa，ua，xb，ub，t)％a表示下邊界，b表示上邊界pa=[0；ua(2)]；qa=[1；0]；pb=[ub(1)一1；O]；qb=[O；1]；(3)初值條件改寫為％初值條件函數(shù)functionu0=pdeic(x)u0=[1；0]；(4)最后編寫主調(diào)函數(shù)ticm=0；x=[00.0050.010.050.10.20.50.70.90.950.990.9951]；T=[00.0050.010.050.10.511.52]；sol=pdepe(m，@pdefun，@pdeic，@pdebc，x，t)；U1=sol(：，：，1)；u2=sol(：，：，2)；figuresurf(x，t，u1)title(’ul(x，t)’)xlabel(’Distancex’)ylabel(’Timet’)figuresurf(x，t，u2)tiffe(’u2(x，t)’)xlabel(’Distancex’)ylabel(’Timet’)結(jié)果圖3．2偏微分方程的pdetool解法3.2.1偏微分方程工具箱的功能偏微分方程工具箱(PDEToolbox)提供了研究和求解空間二維偏微分方程問題的一個強大而又靈活實用的環(huán)境。PDEToolbox的功能包括：(1)設(shè)置PDE(偏微分方程)定解問題，即設(shè)置二維定解區(qū)域、邊界條件以及方程的形式和系數(shù)；(2)用有限元法(FEM)求解PDE數(shù)值解；(3)解的可視化。無論是高級研究人員還是初學者，在使用PDEToo1box時都會感到非常方便。只要PDE定解問題的提法正確，那么，啟動MATLAB后，在MATLAB工作空間的命令行中鍵人pdetool，系統(tǒng)立即產(chǎn)生偏微分方程工具箱(PDEToolbox)的圖形用戶界面(GraphicalUserInterface，簡記為GUI)，即PDE解的圖形環(huán)境，這時就可以在它上面畫出定解區(qū)域、設(shè)置方程和邊界條件、作網(wǎng)格剖分、求解、作圖等工作，詳見1．4節(jié)中的例子。我們將在第二章詳細介紹GUI的使用，在第二章給出大量典型例子和應(yīng)用實例。除了用GUI求解PDE外，也可以用M文件的編程計算更為復雜的問題。3.2.2PDEToolbox求解的問題及其背景方程類型PDEToolbox求解的基本方程有橢圓型方程、拋物型方程、雙曲型方程、特征值方程、橢圓型方程組以及非線性橢圓型方程。橢圓型方程：，橢圓型方程：其中是平面有界區(qū)域，c,a,f以及未知數(shù)u是定義在上的實（或復）函數(shù)。拋物型方程：雙曲型方程：.特征值方程：其中d是定義在上的復函數(shù)，是待求特征值。在拋物型方程和雙曲型方程中，系數(shù)c,a,f和d可以依賴于時間t?？梢郧蠼夥蔷€性橢圓型方程：其中c,a,f可以是未知函數(shù)u的函數(shù)。還可以求解如下PDE方程組；利用命令行可以求解高階方程組。對于橢圓型方程，可以用自適應(yīng)網(wǎng)格算法，還能與非線性解結(jié)合起來使用。另外，對于Poission方程還有一個矩形網(wǎng)格的快速求解器。邊界條件(1)Dirichlet條件：(2)Neumann條件：其中是的邊界上的單位外法向量，和是定義在上的函數(shù)。對于特征值問題僅限于齊次條件：和。對于非線性情形．系數(shù)和可以依賴于u；對于拋物型方程和雙曲型方程，系數(shù)可以依賴于時間t。對于方程組情形，邊界條件為(1)Dirichlet條件：(2)Neumann條件：(3)混合邊界條件為:其中的計算要使得Dirichlet條件滿足。在有限元法中，Dirichlet條件也稱為本質(zhì)邊界條件，Neumann條件稱為自然邊界條件。3.3如何使用FDEToolbox3.3.1定解問題的設(shè)置原簡單的辦法是在PDETool上直接使用圖形用戶界面(GUl)。設(shè)置定解問題包括三個步驟：(1)Draw模式：使用CSG(幾何結(jié)構(gòu)實體模型)對話框畫幾何區(qū)域，包括矩形、圓、橢圓和多邊形，也可以將它們組合使用。(2)Boundary模式：在各個邊界段上給出邊界條件，(3)PDE模式：確定方程的類型、系數(shù)c，a，f和dc。也能夠在不同子區(qū)域上設(shè)置不同的系數(shù)(反映材料的性質(zhì))。3.3.2解PDE問題用GUI解PDE問題主要經(jīng)過下面兩個過程（模式）(1)Mesh模式；生成網(wǎng)格．自動控制網(wǎng)格參數(shù)。(2)Solve模式：對于橢圓型方程還能求非線性和自適應(yīng)解。對于拋物型和雙曲型力程．設(shè)置初始邊值條件后能求出給定t時刻的解。對于特征值問題，能求出給定區(qū)間內(nèi)的特征值；求解后可以加密網(wǎng)格再求解。3.3.3使用Toolbox求解非標準的問題對于非標準的問題?？梢杂肞DEToo1box的函數(shù)?；蛘哂肍EM(有限元法)求解更為復雜的問題。3.3.4計算結(jié)果的可視化從GUI能夠使用Plot模式實現(xiàn)可視化?？梢允褂肅olor,Height和Vector等作圖。對于拋物型和雙曲型方程，還可以生成解的動畫。這些操作通過命令行都很容易實現(xiàn)。3.3.5應(yīng)用領(lǐng)域在應(yīng)用界面提供了丁如下應(yīng)用領(lǐng)域．結(jié)構(gòu)力學——平面應(yīng)力問題．結(jié)構(gòu)力學——平面應(yīng)變問題．靜電場問題．靜磁場問題．交流電磁場問題．直流導體介質(zhì)問題．熱傳導問題3.4解偏微分方程的一個例子解Poisson方程，邊界條件為齊次Dirichlet類型。第一步：啟動MATLABl,鍵入pdetool，按回車鍵確定便可啟動GUI，然后在Options菜單下選擇Grid命令，打開柵格,柵格的使用，能使用戶容易確定所繪圖形的大小，如圖1—11--1第二步：分步完成平面幾何造型：R1-C1-E1+R2+C2。用菜單或快捷工具，分別畫矩形R1、矩形R2、橢圓E1、圓C1、圓C2。畫圓時，首先選中橢圓工具，按鼠標右鍵并拖動即可、或者在按ctrI的同時，拖動鼠標也可繪制圓。然后在Setformula欄，進行編輯并用算術(shù)運算將將圖形對象名稱連接起來，刪除默認的表達式鍵入R1-C1-E1+R2+C2，按等號健得到所需圖形。若需要，還可進行儲存．形成M文件。選擇Boundary菜單中BoundaryMode命令，進入邊界模式。單擊Boundary菜單中RemoveA11SubdomainBorders選項，去除子域邊界。如果想將幾何信息和邊界信息進行存儲，應(yīng)選擇Boundary菜單中的ExPortDecomposedGeometry．BoundaryCond’s…命令，將它們分別儲存于g,b變量中,通過MATLAB形成M文件。第三步：選取邊界．單擊Boundary菜單中SpecifyBounddyConditions…選項，打開Boundaryconditlons對話框，輸入邊界條件，如圖1—4。本例取缺省條件。即將全部邊界設(shè)為齊次Dirichlet條件，邊界顏色顯示為紅色。第四步：選擇PDE菜單中PDEMode命令，進入PDE模式。單擊PDE菜單中PDESpecification…選項，打開PDE對話框，設(shè)置方程類型。本例取缺省設(shè)置，類型為橢圓型，參數(shù)c,a,f分別為1，0，10。第五步：選擇Mesh菜單中InitializeMesh命令，進行網(wǎng)格剖分。第六步：選擇Mesh菜單中RefineMesh命令，對網(wǎng)格加密。第七步：選擇Solve菜單中So1vePDE命令，解偏微分方程并顯示圖形解。第八步：單擊Plot菜單中Parameters…選項，打開Plotselection對話框，選中Color,Height(3—DPlot)和Showmesh三項。然后單擊Plot按鈕，顯示三維圖形解。第九步：如果要畫等值線圖和矢量場圖，單擊Plot菜單中Parameters…選項，打開PlotSelection對話框．選中Contour和Arrows兩項。然后單擊P1ot按鈕，可顯示解的等值線圖和矢量場圖。3.4.1PDEToolbox菜單1.File菜單（如圖1-1）圖1-1New新建一個幾何結(jié)構(gòu)實體模型（ConstructiveSolidGeomery,簡記為CSG），默認文件名為“Untitled”。Open…從硬盤裝載M文件Save將在GUI內(nèi)完成的成果儲存到一個M文件中。SaveAs…將在GUI內(nèi)完成的成果儲存到另外一個M文件中。Print…將PDE工具箱完成的圖形送到打印機內(nèi)進行硬拷貝。Exit退出PDE工具圖形用戶界面。2、Edit菜單（如圖1-2）圖1-2Undo在繪制多邊形時退回到上一步操作。Cut將已選實體剪切到剪貼板上。Copy將已選實體拷貝到剪貼板上。Paste…將剪貼板上的實體粘貼到當前幾何結(jié)構(gòu)實體模型中。Clear刪除已選的實體。SelectAll選擇當前幾何結(jié)構(gòu)實體造型CSG中的所有實體及其邊界和字域。3Options菜單（如圖1-3）圖1-3Grid繪圖時打開或關(guān)閉柵格。GridSpacing…調(diào)整柵格的大小。Snap打開或關(guān)閉捕捉柵格功能。AxesLimits…設(shè)置繪圖軸的坐標范圍。AxesEqual打開或關(guān)閉繪圖方軸。TurnoffToolbarHelp關(guān)閉工具欄按鈕的幫助信息。Zoom打開或關(guān)閉圖形縮放功能。Application選擇應(yīng)用的模式。Refresh重新顯示PDE工具箱中的圖形實體。4、Draw菜單（如圖1-4）圖1-4DrawMode進入繪圖模式。Rectangle/square以角點方式畫矩形/方行（Ctrl+鼠標）。Rectangle/square（centered）以中心方式畫矩形/方行（Ctrl+鼠標）。Ellipse/circle以矩陣角點方式畫橢圓/圓（Ctrl+鼠標）。Ellipse/circle（centered）以中心方式畫橢圓/圓（Ctrl+鼠標）。Polygon畫多邊形，單擊鼠標右鍵可封閉多邊形。Rotate…旋轉(zhuǎn)已選的圖形。ExportGeometryDescription,SetFormula，Labels…將幾何描述矩陣gd、公式設(shè)置字符sf和標識空間矩陣ns輸出到主工作空間去。單擊Draw菜單中Rotate…選項，可打開Rotate比對活框，通過輸入旋轉(zhuǎn)的角度，可使選擇的物體按輸入的角度逆時針旋轉(zhuǎn)。旋轉(zhuǎn)中心的選擇如果缺省，則為圖形的質(zhì)心，也可以輸入旋轉(zhuǎn)中心坐標。5、Boundary菜單（如圖1-5）圖1-5BoundaryMode進入邊界模式。SpecifyBoundaryConditions…對于已選的邊界輸入條件，如果沒有選擇邊界，則邊界條件適用于所有的邊界。ShowEdgeLabels顯示邊界區(qū)域標識開關(guān)，其數(shù)據(jù)是分解幾何矩陣的列數(shù)。ShowSubdomainLabels顯示子區(qū)域標識開關(guān)，其數(shù)據(jù)是分解幾何矩陣中的子域數(shù)值。RemoveSubdomainBorder當圖形進行布爾運算時，刪除已選取的子域邊界。RemoveAllSubdomainBorders當圖形進行布爾運算時，刪除所有的子域邊界。ExportDecomposedGeometry，BoundaryCond’s…將分解幾何矩陣g、邊界條件矩陣b輸出到主工作空間。選擇Boundary菜單中SpecifyBoundaryConditions．命令可定義邊界條件。在打開的Boundarycondition對話框，可對已選的邊界輸入邊界條件。共有如下三種不同的條件類型：NeMmann條件這里邊界條件是由方程系數(shù)q和g確定的，在方程組的情況下(換成方程組模式)，q是2ⅹ2矩陣，g是2x1矢量。Dirichlet條件u定義在邊界上，邊界條件方程是價h*u=r，這里h是可以選樣的權(quán)因子(通常為1)。在方程組情況下，h是2x2矩陣，r是2xl矢量，混合邊界條件(僅適合于方程組情形)它是Dirichlet和Neumann的混合邊界條件，q是2x2矩陣，g是2x1矢量，h是1x2矢量，r是一個標量。6、PDE菜單（如圖1-6）圖1-6PDEMode進入偏微分方程模式。ShowSubdomainLabels顯示子區(qū)域標識開關(guān)。PDESpecification…調(diào)整PDE參數(shù)和類型。ExportPDECoefficients…將當前PDE參數(shù)c,a,f,d輸出到主工作空間，其參數(shù)變量為字符類型。7、Mesh菜單（如圖1-7）圖1-7MeshMode輸入網(wǎng)格模式。InitializeMesh建立和顯示初始化三角形網(wǎng)格。RefineMesh加密當前三角型網(wǎng)格。JiggleMesh優(yōu)化網(wǎng)格。UndoMeshChange退回上一次網(wǎng)格操作。DisplayTriangleQuality用0~1之間數(shù)字化的顏色顯示三角形網(wǎng)格的質(zhì)量，大于0.6的網(wǎng)格可接受的。ShowNodeLabels顯示網(wǎng)格節(jié)點標識開關(guān)，節(jié)點標識數(shù)據(jù)是點矩陣p的列。ShowTriangleLabels顯示三角形網(wǎng)格標識開關(guān)，三角形網(wǎng)格標識數(shù)據(jù)是三角形矩陣t的列。Parameters…修改網(wǎng)格生成參數(shù)。ExportMesh輸出節(jié)點矩陣p、邊界矩陣e和三角形矩陣t到主工作空間。8、Solve菜單（如圖1-8）圖1-8SolvePDE對當前的幾何結(jié)構(gòu)實體CSG、三角形網(wǎng)格和圖形解偏微分方程。Parameters…調(diào)整PDE的參數(shù)。ExportSolution…輸出PDE的解矢量u。如果可行，將計算的特征值1輸出到主工作空間。9、Plot菜單（如圖1-9）圖1-9PlotSolution顯示圖形解。Parameters…打開繪圖方式對話框。ExportMovie…如果動畫被錄制了，則動畫矩陣M將輸出到主工作空間。10、Window菜單從Window菜單項，可選擇當前打開的所有的MATLAB圖形窗口，被選擇的窗臺移至前臺。11、Help菜單Help…顯示幫助信息About…顯示版本信息主菜單下是工具欄，工具欄中喊有許多工具圖標按鈕，可提供快速、便捷的操作方式。從左到右5個按鈕為繪圖模式按鈕，緊接著的6個為邊界、網(wǎng)格、解方程和圖形顯示控制功能按鈕，最右邊的為圖形縮放功能鍵。（如圖1-10）圖1-10以角點方式畫矩形/方行（Ctrl+鼠標）。以中心方式畫矩形/方行（Ctrl+鼠標）。以矩形角點長軸方式畫橢圓/圓（Ctrl+鼠標）。以中心方式畫橢圓/圓（Ctrl+鼠標）。畫多邊形，按右鍵可封閉多邊形。進入邊界模式。打開PDESpecification（偏微分方程類型）對話框。初始化三角形網(wǎng)格。加密三角形網(wǎng)格。解偏微分方程。打開PlotSelection對話框，確定后給出解的三維圖形。為顯示縮放切換按鈕。求解PDE問題主要有兩種方法，一種是使用圖形用戶界面，另一種是采用命令行編程。前者直觀簡便，而后者更為靈活。3.4.2求解橢圓方程及拋物型方程的例子例1：單位圓上的Poisson方程邊值問題：這一問題的精確解為：若使用圖形用戶界面（GraphicalUserInterface,簡記為GUI），則首先在MATLAB的工作窗口中鍵入pdetool，按回車鍵確定，于是出現(xiàn)PDEToolbox窗口。如果需要坐標網(wǎng)格，單擊Options菜單下的Grid選項即可。下面分步進行操作。（i）畫區(qū)域圓單擊工具，大致在（0，0）位置單擊鼠標右鍵同時拖拉鼠標到適當位置松開，繪制圓。為了保證所繪制的圓是標準的單位圓，在所繪圓上雙擊，打開ObjectDialog對話框，精確地輸入圓心坐標X-center為0、Y-cebter為0及半徑Radius為1，然后單擊OK按鈕，這樣單位遠已畫好。（ii）設(shè)置邊界條件單擊工具，圖形邊界變紅，逐段雙擊邊界，打開BoundaryCondition對話框，輸入邊界條件。對于同一類型的邊界，可以按Shift鍵，將多個邊界同時選擇，統(tǒng)一設(shè)置邊界條件。本題選擇Dirichlet條件，輸入h為1，r為0，然后單擊OK按鈕。也可以單擊Boundary菜單中SpecifyBoundaryConditions…選項，打開BoundaryCondition對話框輸入邊界條件，如圖2-1。（iii）設(shè)置方程單擊PDE菜單中PDESpecification…選項，打開PDESpecification對話框，選項方程類型。本題單擊Elliptic，輸入c為1，a為0，f為1，然后單擊OK按鈕，如圖2-2。圖2-1圖2-2（iv）網(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ù)變?yōu)?41個，三角形單元數(shù)為1016個，如此還可繼續(xù)加密。（v）解方程單擊工具，或者單擊Solve菜單中Solve菜單中SolvePDE選項，可顯示方程色彩解。如果單擊Plot菜單中Parameters…選項，出現(xiàn)PlotSelection對話框，如圖2-3，從中可以選擇Color，Contour，Arrows,Deformedmesh，Height（3-Dpolt），還可以設(shè)置等值線的數(shù)目等。本例中選擇Color，Contour，Height（3-Dpolt）和Showmesh四項，然后單擊Plot按鈕，方程的圖形解如圖2-4所示。除了作定解問題解u的圖形外，也可以作|gradu|，|cgradu|等圖形。圖2-3圖2-4（vi）與精確解作比較單擊Plot菜單中Parameters…選項，打開PlotSelection對話框，在Height（3-Dplot）行Property下拉框中選userentry，且在該行的Userentry輸入框中鍵入u-（1-x.^2-y.^2）/4,單擊Ｐl(wèi)ot按鈕就可以看到解的絕對誤差圖形，如圖２-５．可見在邊界處誤差為０。圖2-5（vii）輸出網(wǎng)格節(jié)點的編號、單元編號以及節(jié)點坐標單擊Mesh菜單中ShowNodeLabels選項，再單擊工具或，即可顯示節(jié)點編號。若要輸出節(jié)點坐標，只需單擊Mesh菜單中ExportMesh…選項，這時打開的Export對話框中默認值為p，e，t，這里p，e，t分別表示points（點）、edges（邊）、triangles（三角形）數(shù)據(jù)的變量，單擊OK按鈕。然后在MATLAB命令窗口鍵入p，按回車鍵確定，即可顯示出節(jié)點按編號排列的坐標（二維數(shù)組）；鍵入e，按回車鍵，則顯示邊界線段數(shù)據(jù)矩陣（7維數(shù)組）；輸入t，按回車鍵，則顯示三角形單元數(shù)據(jù)矩陣（4維數(shù)組）。（viii）輸出解的數(shù)值單擊Solve菜單中ExportSolution…選項，在打開的Export對話框中輸入u，單擊OK按鈕確定。再在MATLAB命令窗口中輸入u，按回車確定，即顯示按節(jié)點編號排列的解的數(shù)值。我們也可以用MATLAB程序求解PDE問題，同時顯示解的圖形：[p,e,t]=initmesh(‘circleg’,’hmax’,1);Error=[];err=1;Whileerr>0.001,[p,e,t]=refinemesh(‘circleg’,p,e,t);U=assempde(‘circleb1’,p,e,t,,1,0,1);Exact=(1-p(1,:).^2-p(2,:).^2)’/4;Err=norm(u-exact,’inf’);Error=[errorerr];EndPdemesh(p,e,t)Pdesurf(p,t,u)Pdesurf(p,t,u-exact)通過命令行鍵入help+命令函數(shù)，如helppdemesh，按回車鍵，可以調(diào)入有關(guān)命令函數(shù)的定義、參數(shù)格式等幫助信息。例2：求解拋物型方程的例子考慮一個帶有矩形孔的金屬板上的熱傳導問題。板的左邊保持在100，板的右邊熱量從板向環(huán)境空氣定常流動，其他邊及內(nèi)孔邊界保持絕緣。初始時板的溫度為0，于是概括為如下定解問題：域的外邊界頂點坐標為（-0.5，-0.8），（0.5，-0.8），（0.5，0.8），（-0.5，0.8）。內(nèi)邊界頂點坐標為（-0.005，-0.4），（0.05，-0.4），（0.05，0.4），（-0.05，0.4）。使用GUI求解這一問題。在PDEToolbox窗口的工具欄中選擇GenericScalar模式。（i）區(qū)域設(shè)置單擊工具，在窗口拖拉出一個矩形，雙擊矩形區(qū)域，在ObjectDialog對話框中輸入Left為-0.5，Bottom為-0.8，Width為1，Height為1.6，單擊OK按鈕，顯示矩形區(qū)域R1。用同樣方法作內(nèi)孔R2，只要設(shè)置Left=-0.05，Bottom=-0.4，Width=0.1，Height=0.8即可。然后在Setformula欄中鍵入R1-R2。（ii）設(shè)置邊界條件單擊，使邊界變紅色，然后分別雙擊每段邊界，打開BoundaryCvondition對話框，設(shè)置邊界條件。在左邊界條件。在左邊界上，選擇Dirichlet條件，輸入h為1，r為100；右邊界上，選擇Neumann條件，輸入g為-1，q為0；其他邊界上，選擇Neumann條件，輸入g為0，q為0。（iii）設(shè)置方程類型單擊，打開PDESpecification對話框，設(shè)置方程類型為Parabolic（拋物型），d=1，c=1，a=0，f=0，單擊OK按鈕。（iv）網(wǎng)格剖分單擊，或者加密網(wǎng)格，單擊。（v）初值和誤差的設(shè)置單擊Solve菜單中Parameters…選項，打開SolveParameters對話框，輸入Time為0:5，u（t0）為0，Relativetolerance為0.01，Absolutetolerance為0.001，然后單擊OK按鈕。（vi）數(shù)值解的輸