偏微分方程的MATLAB解法

李 明

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

1 前言

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

2 Matlab解偏微分方程

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

2.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：是一個(gè)三維數(shù)組，sol（：，：，i）表示ui的解，換句話說uk對(duì)應(yīng)x（i）和t（j）時(shí)的解為sol（i，j，k），通過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)注形式，則原方程可以改寫為

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

（2）邊界條件改寫為

%邊界條件函數(shù)

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

（3）初值條件改寫為

%初值條件函數(shù)

（4）最后編寫主調(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)系中畫出偏微分方程的大致定解區(qū)域。

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

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

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

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

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

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

（10）如要畫出數(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命令，打開柵格。柵格的使用，能使用戶容易確定所繪圖形的大小。

第二步：分步完成平面幾何造型：R1-C1-E1+R2+C2。用菜單或快捷工具，分別畫矩形R1、矩形R2、橢圓E1、圓C1、圓C2．畫圓時(shí)，首先選中橢圓工具，按鼠標(biāo)右鍵并拖動(dòng)即可，或者在按Ctrl的同時(shí)，拖動(dòng)鼠標(biāo)也可繪制圓。然后在Set formula欄，進(jìn)行編輯并用算術(shù)運(yùn)算符將圖形對(duì)象名稱連接起來，或刪除默認(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變量中，通過MATLAB形成M文件。

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

第四步：選擇PDE菜單中PDE Mode命令，進(jìn)入PDE模式。單擊PDE菜單中PDE Specification…選項(xiàng)，打開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)，打開Plot Selection 對(duì)話框，選中 Color，Height（3-D Plot）和 Show mesh三項(xiàng)．然后單擊Plot按鈕，顯示三維圖形解。

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

3 結(jié)語

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

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

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

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