4.8.2? 偏微分方程

在自然科學的很多領域內(nèi)，都會遇到微分方程初值問題，特別是偏微分方程，它的定解問題是描述自然界及科學現(xiàn)象的最重要的工具。可以說，幾乎自然界和各種現(xiàn)象都可以通過微分方程(特別是偏微分方程)來描述。

MATLAB提供了一個專門用于求解偏微分方程的工具箱PDE Toolbox。本小節(jié)僅介紹一些最簡單、經(jīng)典的偏微分方程，如橢圓型、雙曲型、拋物型等偏微分方程，并給出求解方法。用戶可以從中了解其解題的基本方法，從而解決類似的問題。

1．橢圓型問題

assempde函數(shù)是PDE工具箱中的一個基本函數(shù)，它使用有限元法組合PDE問題。該函數(shù)用來有選擇地生成PDE問題的解。可以用assempde函數(shù)求解下面的標量橢圓型問題：

或系統(tǒng)橢圓型問題：

對于標量的情況，用解的列向量代表解矢量u，列矢量中的值對應于p的對應節(jié)點處的解。對于具有np個節(jié)點的N維系統(tǒng)，u1的前np行描述u的第1個元素，接下來的np行描述u的第2個元素，依次類推。這樣，u的元素就作為N塊節(jié)點行放到u中。assempde函數(shù)的調(diào)用語法如下。

(1)u=assempde(b,p,e,t,c,a,f)： 通過在線性方程組中剔除Dirichlet邊界條件來組合和求解PDE問題。

(2)[K,F]=assempde(b,p,e,t,c,a,f)：?通過剛性位移近似Dirichlet邊界條件來組合和求解PDE問題。K和F分別為剛性矩陣和右邊項。PDE問題的有限元解為u1=K\F。

(3)[K,F,B,ud]=assempde(b,p,e,t,c,a,f)：通過從線性方程組剔除Dirichlet邊界條件來組合PDE問題。u1=K\F，則返回非Dirichlet點上的解。完整的PDE問題可以通過MATLAB中的表達式u=B*u1+ud求解。

(4)[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f)：給出PDE問題的分離表示。

(5)u=assempde(K,M,F,Q,G,H,R)：將PDE問題的分離表示轉(zhuǎn)換為單一矩陣或矢量的形式，然后通過從線性方程組中剔除Dirichlet邊界條件來求解。

(6)[K1,F1]=assempde(K,M,F,Q,G,H,R)：用很大的常數(shù)修改Dirichlet邊界條件，從而將PDE問題的分離表示轉(zhuǎn)換為單一矩陣或矢量的形式。

(7)[K1,F1,B,ud]=assempde(K,M,F,Q,G,H,R)：從線性方程組中剔除Dirichlet邊界條件，從而將PDE問題的分離表示轉(zhuǎn)換為單一矩陣或矢量的形式。

b描述PDE問題的邊界條件。b也可以是邊界條件矩陣或邊界M文件的文件名。PDE問題幾何模型由網(wǎng)絡數(shù)據(jù)p,e,t描述。網(wǎng)格數(shù)據(jù)的生成可以查詢help文檔中的initmesh函數(shù)。

系數(shù)c，a，d，f可以通過多種方式給定。這些系數(shù)也可以與時間t相關，在assempde函數(shù)中可以看到所有選項的列表。

【例4-47】? 求解L型薄膜的方程，為Dirichlet邊界條件。最后繪圖顯示結(jié)果。

>> [p,e,t]=initmesh('lshapeg','hmax',0.2); ?? %生成初始三角形網(wǎng)格，hmax指網(wǎng)格大小

>> [p,e,t]=refinemesh('lshapeg',p,e,t);???? ? %? 將初始的三角形網(wǎng)格細化

>> u=assempde('lshapeb',p,e,t,1,0,1);?????? ? %? 求解方程

>> pdesurf(p,t,u)?????????????????????????????? ?%? 繪制結(jié)果圖形

lshapeg和lshapeb分別為表示對象幾何模型和邊界條件的M文件，為工具箱自帶文件。initmesh函數(shù)和 refinemesh函數(shù)分別對網(wǎng)格模型進行初始化和細化。pdesurf函數(shù)繪制解的表面圖。L型薄膜泊松方程的解如圖4-15所示。

2．拋物型問題

MATLAB提供了parabolic函數(shù)求解標量拋物型問題：

圖4-15? L型薄膜泊松方程的解

或系統(tǒng)PDE問題：

該函數(shù)的調(diào)用語法如下。

(1)u1=parabolic(u0,tlist,g,b,p,e,t,c,a,f,d)：p,e,t為網(wǎng)格數(shù)據(jù)，b為邊界條件，初值為u0。

(2)u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d,rtol,atol)：atol和rtol為絕對和相對容限。

(3)u1=parabolic(u0,tlist,K,F,B,ud,M)：生成下面PDE問題的解。

的初始值為u0。

【例4-48】? 求解熱傳導方程：。其中。在的區(qū)域內(nèi)令，在其他區(qū)域內(nèi)令。使用Dirichlet邊界條件。要求計算時間linspace(0,0.1,20)處的解。

%? 生成初始三角形網(wǎng)格數(shù)據(jù)

>> [p,e,t]=initmesh('squareg');????????????? ? %? 參數(shù)squareg指計算區(qū)域是方形的

>> [p,e,t]=refinemesh('squareg',p,e,t);????? ?%? 將初始的三角形網(wǎng)格細化

>> u0=zeros(size(p,2),1);

>> ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);? ? %? 查找區(qū)域內(nèi)的元素

>> u0(ix)=ones(size(ix));????????????????????? ?%? 令

>> tlist=linspace(0,0.1,20);????????????????? ?%? 時間列表

>> u1=parabolic(u0,tlist,'squareb1',p,e,t,1,0,1,1);?? %? 求解偏微分方程

本例首先調(diào)用initmesh函數(shù)生成偏微分方程的初始網(wǎng)格，然后調(diào)用parabolic函數(shù)求解偏微分方程，運行結(jié)束將顯示如下信息：

96 successful steps

0 failed attempts

194 function evaluations

1 partial derivatives

20 LU decompositions

193 solutions of linear systems

具體的結(jié)果u1是一個665*20的矩陣，這里就略去不顯示了。

3．雙曲型問題

MATLAB提供了hyperbolic函數(shù)來求解標量雙曲型問題：

或系統(tǒng)雙曲型問題：

hyperbolic函數(shù)的調(diào)用語法如下。

u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d,rtol,atol)：p,e,t為網(wǎng)格數(shù)據(jù)，b為邊界條件，u0為初值，初始導數(shù)為ut0。atol和rtol為絕對和相對容限。

u1=hyperbolic(u0,ut0,tlist,K,F,B,ud,M)生成下面PDE問題的解。

u的初始值為u0和ut0。

在MATLAB命令窗口輸入：

>> [p,e,t]=initmesh('squareg');????????????? ????????????? %?? 生成初始三角形網(wǎng)格

>> x=p(1,:)';

>> y=p(2,:)';

>> u0=atan(cos(pi/2*x));

>> ut0=3*sin(pi*x).*exp(cos(pi*y));

>> tlist=linspace(0,5,31);??????????????????? ????????????? %? 時間列表

>> uu=hyperbolic(u0,ut0,tlist,'squareb3',p,e,t,1,0,0,1);? ?%? 求解方程

本例首先調(diào)用initmesh函數(shù)生成偏微分方程的初始網(wǎng)格，然后調(diào)用hyperbolic函數(shù)求解偏微分方程，運行結(jié)束將顯示如下信息：

462 successful steps

70 failed attempts

1066 function evaluations

1 partial derivatives

156 LU decompositions

1065 solutions of linear systems

具體的結(jié)果uu是一個177*31的矩陣，這里就略去不顯示了。

總結(jié)

以上是生活随笔為你收集整理的matlab画微分方程的矢量场图_MATLAB偏微分方程的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。