偏微分方程的MATLAB解法

李 明

（長(zhǎng)江大學(xué) 信息與數(shù)學(xué)學(xué)院，湖北 荊州 434023）

1 前言

偏微分方程定解問(wèn)題有著廣泛的應(yīng)用背景。人們用偏微分方程來(lái)描述、解釋或者預(yù)見(jiàn)各種自然現(xiàn)象，并用于科學(xué)和工程技術(shù)的各個(gè)領(lǐng)域[1]。然而，對(duì)于廣大應(yīng)用工作者來(lái)說(shuō)，從偏微分方程模型出發(fā)，使用有限元法或有限差分法求解都要耗費(fèi)很大的工作量，才能得到數(shù)值解。現(xiàn)在，MATLAB PDE Toolbox已實(shí)現(xiàn)對(duì)于空間二維問(wèn)題高速、準(zhǔn)確的求解過(guò)程。

2 Matlab解偏微分方程

MATLAB提供了兩種方法[2]解決PDE問(wèn)題，一是pdepe函數(shù)[3]，它可以求解一般的PDEs，具有較大的通用性，但只支持命令行形式調(diào)用。二是PDE工具箱[4]，可以求解特殊PDE問(wèn)題，但有較大的局限性，比如只能求解二階PDE問(wèn)題，并且不能解決偏微分方程組。它提供了GUI界面，可以從繁雜的編程中解脫出來(lái)，同時(shí)還可以通過(guò)FileSave As直接生成M代碼。

2.1 函數(shù)解法

pdepe函數(shù)介紹

它的調(diào)用格式為 sol=pdepe（m，@pdefun，@pdeic，@pdebc，x，t）

輸入?yún)?shù)： @pdefun：是PDE的問(wèn)題描述函數(shù)

@pdebc：是PDE的邊界條件描述函數(shù)

@pdeic：是PDE的初值條件

輸出參數(shù)：sol：是一個(gè)三維數(shù)組，sol（：，：，i）表示ui的解，換句話說(shuō)uk對(duì)應(yīng)x（i）和t（j）時(shí)的解為sol（i，j，k），通過(guò)sol，我們可以使用pdeval（）直接計(jì)算某個(gè)點(diǎn)的函數(shù)值。

實(shí)例講解

例：試求解下面的偏微分

其中，F(xiàn)（x）=e5.73x-e-11.46x，且滿足初始條件u1（x，0）＝1，u2（x，0）＝0 及邊界條件

解：（1）對(duì)照給出的偏微分方程，根據(jù)標(biāo)注形式，則原方程可以改寫(xiě)為

%目標(biāo)PDE函數(shù)

（2）邊界條件改寫(xiě)為

%邊界條件函數(shù)

function[pa，qa，pb，qb]=pdebc（xa，ua，xb，ub，t）%a表示下邊界，b表示上邊界

（3）初值條件改寫(xiě)為

%初值條件函數(shù)

（4）最后編寫(xiě)主調(diào)函數(shù)

圖1 結(jié)果圖

2.2 偏微分方程的pdetool解法

2.2.1 pdetool介紹pdetool提供的用戶圖形界面（GUI）解法的使用步驟如下：

（1）在 Matlab命令窗口運(yùn)行 pdetool，出現(xiàn) PDE Toolbox界面。

（2）用鼠標(biāo)點(diǎn)一下工具欄上的“PDE”按鈕，在彈出的對(duì)話框中定義偏微分方程。

（3）用鼠標(biāo)點(diǎn)一下工具欄上的區(qū)域按鈕，在下面的坐標(biāo)系中畫(huà)出偏微分方程的大致定解區(qū)域。

（4）雙擊（3）中畫(huà)出的大致區(qū)域，在彈出的對(duì)話框中精確定位定解區(qū)域。

（5）用鼠標(biāo)點(diǎn)一下工具欄上的邊界按鈕“墜Ω”，畫(huà)出區(qū)域的邊界。

（6）雙擊坐標(biāo)系中的區(qū)域邊界，定義偏微分方程的邊界條件。

（7）用鼠標(biāo)點(diǎn)工具欄上的剖分按鈕，對(duì)求解區(qū)域進(jìn)行剖分。

（8）如果求拋物型或雙曲型方程的數(shù)值解，還需要通過(guò)“solve”菜單下的“parameters…”選項(xiàng)設(shè)置初值條件。

（9）用鼠標(biāo)點(diǎn)一下工具欄上的“＝”按鈕，就畫(huà)出偏微分方程數(shù)值解的圖形。通過(guò)“solve”菜單下的“Export Solution…”選項(xiàng)可以把數(shù)值解輸出到Matlab的工作間。

（10）如要畫(huà)出數(shù)值解的三維圖形，需要設(shè)置“plot”菜單下的“parameters…”選項(xiàng)。

2.2.2 實(shí)例講解

例：解Poisson方程-△u=f，邊界條件為齊次Dirichlet類型。

第一步：?jiǎn)?dòng)MATLAB，鍵入pdetool，按回車鍵確定便可啟動(dòng)GUI，然后在Options菜單下選擇Grid命令，打開(kāi)柵格。柵格的使用，能使用戶容易確定所繪圖形的大小。

第二步：分步完成平面幾何造型：R1-C1-E1+R2+C2。用菜單或快捷工具，分別畫(huà)矩形R1、矩形R2、橢圓E1、圓C1、圓C2．畫(huà)圓時(shí)，首先選中橢圓工具，按鼠標(biāo)右鍵并拖動(dòng)即可，或者在按Ctrl的同時(shí)，拖動(dòng)鼠標(biāo)也可繪制圓。然后在Set formula欄，進(jìn)行編輯并用算術(shù)運(yùn)算符將圖形對(duì)象名稱連接起來(lái)，或刪除默認(rèn)的表達(dá)式直接鍵入R1-C1-E1+R2+C2。若需要，還可進(jìn)行儲(chǔ)存，形成M文件。

選擇Boundary菜單中Boundary Mode命令，進(jìn)入邊界模式．單擊Boundary菜單中Remove All Subdomain Borders選項(xiàng)，去除子域邊界．如果想將幾何信息和邊界信息進(jìn)行存儲(chǔ)，應(yīng)選擇Boundary菜單中的Export Decomposed Geometry，Boundary Cond’s…命令，將它們分別存儲(chǔ)于g，b變量中，通過(guò)MATLAB形成M文件。

第三步：選取邊界，單擊Boundary菜單中Secify Boundary Conditions…選項(xiàng)，打開(kāi)Boundary Conditions對(duì)話框，輸入邊界條件．本例取缺省條件，即將全部邊界設(shè)為齊次Dirichlet條件，邊界顏色顯示為紅色。

第四步：選擇PDE菜單中PDE Mode命令，進(jìn)入PDE模式。單擊PDE菜單中PDE Specification…選項(xiàng)，打開(kāi)PDE Specification對(duì)話框，設(shè)置方程類型。本例取缺省設(shè)置，類型為橢圓型，參數(shù) c，a，f分別為 1，0，10．

第五步：選擇Mesh菜單中Initialize Mesh命令，進(jìn)行網(wǎng)格剖分。

第六步：選擇Mesh菜單中Refine Mesh命令，對(duì)網(wǎng)格加密。

第七步：選擇Solve菜單中Solve PDE命令，解偏微分方程并顯示圖形解。

第八步：?jiǎn)螕鬚lot菜單中Parameters…選項(xiàng)，打開(kāi)Plot Selection 對(duì)話框，選中 Color，Height（3-D Plot）和 Show mesh三項(xiàng)．然后單擊Plot按鈕，顯示三維圖形解。

第九步：如果要畫(huà)等值線圖和矢量場(chǎng)圖，單擊Plot菜單中Parameters…選項(xiàng)，打開(kāi)Plot Selection對(duì)話框，選中Contour和Arrows兩項(xiàng)。然后單擊Plot按鈕，可顯示解的等值線圖和矢量場(chǎng)圖。

3 結(jié)語(yǔ)

偏微分方程是一門(mén)實(shí)用性很強(qiáng)的學(xué)科。對(duì)于理論研究和實(shí)際應(yīng)用中提出的偏微分方程問(wèn)題，由于其邊界和邊界條件復(fù)雜等原因，尋求解析解是非常困難甚至不可能的，利用計(jì)算機(jī)研究相應(yīng)問(wèn)題的數(shù)值解法是十分必要的。編程實(shí)現(xiàn)從偏微分方程數(shù)值求解全過(guò)程需要很好的理論基礎(chǔ)和編程技巧，而偏微分方程工具箱提供了研究和求解空間二維偏微分方程問(wèn)題的一個(gè)強(qiáng)大而又靈活實(shí)用的環(huán)境.借助于這個(gè)工具，我們可以從繁瑣，共性的求解步驟中解脫出來(lái)而專注于問(wèn)題的核心即問(wèn)題的描述，定義及簡(jiǎn)化，邊界條件的確定，求解方法和精度控制的選擇等，大大提高了求解效率。

[1]陸君安，尚濤，謝進(jìn)，等.偏微分方程的 MATLAB 解法[M].武漢：武漢大學(xué)出版社，2001.

[2]王定江.應(yīng)用偏微分方程[M].杭州：浙江大學(xué)出版社，2007.

[3]王家文，王皓，劉海.MATLAB7.0編程基礎(chǔ)[M].北京：機(jī)械工業(yè)出版社，2005.