1、MATLAB求解PDE問題（1）概述、例子已有2363次閱讀2010-10-12 14:57 |個(gè)人分類:生活點(diǎn)滴|系統(tǒng)分類:科研筆記 |關(guān)鍵詞:MATLAB PDE Toolbo橢圓型方程 有限元方法MATLAB PDEofolbox提供利用有限元方法求解偏微分方程的 GUI以及相應(yīng) 的命令行函數(shù)。利用該工具箱可以求解橢圓型方程、拋物型方程、雙曲型方程、 特征值方程以及非線性方程。PDE Toolbox的功能非常強(qiáng)大，網(wǎng)上有許多利用 PDE Toolbox解決各種物理問題的論文，還有專門介紹工具箱的參考書。網(wǎng)上的例子雖然很多，但是大部分是介紹 PDE工具箱自帶的一些例子，這些 例子中解的區(qū)域

2、，邊界條件是PDE工具箱已經(jīng)編寫好的，直接調(diào)用就可以。對(duì)于 該如何自己設(shè)定求解區(qū)域及邊界條件， 卻很少有人涉及。網(wǎng)上搜索發(fā)現(xiàn)只有 劉平 在博客中詳細(xì)介紹過求解區(qū)域的設(shè)定。下面以一個(gè)橢圓型方程的例子來詳細(xì)說明 求解的各個(gè)步驟，希望對(duì)大家能有所幫助。設(shè)要求如下形式的橢圓方程的解：宀幾+ 2(."計(jì)-4 = 0邊界條件為y=o./= r.2 1仝CA'按照PDE的要求，將方程化為標(biāo)準(zhǔn)形式-V -（其中c- 1. 7 = 0 * f -求解后的圖像如下，第一幅圖是解的圖像，第二幅是計(jì)算誤差。從第二幅圖 可以看到，計(jì)算的最大誤差是10-3方量級(jí)。通過這個(gè)例子我們可以基本掌握 PDE求解

3、偏微分方程的步驟和方法，后面 我將詳細(xì)介紹如何設(shè)置區(qū)域及邊界條件。掌握了區(qū)域和邊界條件的設(shè)定，就可以 輕松求解遇到的偏微分方程了。圖后是附帶的matlab命令以及注釋，并提供 m文件附件下載，下載后解壓即可。希望能對(duì)大家有所幫助。x IDx ID下面是編寫的求解上述方程的 matlab語(yǔ)句及說明：g='mygeom' b='mybo un d'定義區(qū)域，邊界條件。mygeom是定義區(qū)域的子函數(shù)名，函數(shù)名可根據(jù)自己的需 要取定，區(qū)域的確定規(guī)則由pdegeom函數(shù)說明，注意pdegeom函數(shù)只是說明如 何定義區(qū)域，它并不直接確定區(qū)域；mybound是定義邊界條件的子

4、函數(shù)名，與區(qū) 域類似，邊界的確定規(guī)則由函數(shù) pdebou nd確定。后面我會(huì)詳細(xì)介紹區(qū)域和邊界 的取法。p,e,t = in itmesh(g);網(wǎng)格初始化，此處也可以寫成p,e,t = initmesh('mygeom');這樣可以省略上面的語(yǔ) 句p,e,t = refin emesh(g,p,e,t);p,e,t = refin emesh(g,p,e,t);加密網(wǎng)格兩次，需要加密幾次重復(fù)幾次即可，根據(jù)具體問題確定加密次數(shù)U= assempde(b,p,e,t,1,0,'2*(x+y)-4');調(diào)用assempde函數(shù)計(jì)算方程的數(shù)值解，assempde函數(shù)的

5、詳細(xì)用法可以參考MATH 網(wǎng)站或者PDE的使用指南。常用的用法是u,res=assempde(b,p,e,t,c,a,，其中b 為邊界條件，此處也可以寫為mybound', p,e,t,為網(wǎng)格參數(shù)，c,a,f,為方程的參 數(shù)，后面也可以加猜測(cè)值以及各種屬性。pdesurf(p,t,U)grid on;xlabel('x');ylabel('y');zlabel('u')colorbarview(60 30)畫出解的圖形。注意，為了讓結(jié)果更直觀一些，使用view函數(shù)調(diào)整了視點(diǎn)位置。 大家可以自行調(diào)整視角，滿意即可。exact=p(1,:)A

6、2+p(2,:)A2-p(1,:).*p(2,:).*(p(1,:)+p(2,:);exact=exact'figurepdesurf(p,t,U-exact)grid onxlabel('x');ylabel('y');zlabel('error')colorbarview(60 30)由于方程有解析解，我們可以比較數(shù)值計(jì)算的誤差。如果能求得解析解，我們也 不會(huì)設(shè)計(jì)各種方法求數(shù)值解了，因此，這一步在大多數(shù)情況下是用不上的， 這里 只是為了比較計(jì)算結(jié)果，驗(yàn)證計(jì)算的精度。m文件下載MATLAB求解PDE問題(2)確定幾何區(qū)域已有1327次閱

7、讀2010-10-12 18:58 |個(gè)人分類:生活點(diǎn)滴|系統(tǒng)分類:科研筆記 |關(guān)鍵詞:MATLAB PD幾何區(qū)域pdegeom前一篇介紹了如何利用Matlab求解橢圓型方程，下面介紹如何確定求解的 幾何區(qū)域。PDE Toolbox中規(guī)定幾何區(qū)域的 m文件是pdegeom.m。但是pdegeom 并不是一個(gè)可以調(diào)用的函數(shù)，它只是規(guī)定了應(yīng)該何如定義區(qū)域，具體的區(qū)域則要 根據(jù)研究的問題來決定。函數(shù)pdegeom釋義如下：參數(shù)為0個(gè)時(shí)，即沒有參數(shù)時(shí)，返回邊界的段數(shù)；參數(shù)為1個(gè)時(shí)，即只有bs，返回輸出區(qū)域邊界的參變量范圍矩陣 d； 參數(shù)為2個(gè)時(shí)，返回每段邊界長(zhǎng)度為s時(shí)的坐標(biāo)。函數(shù)參數(shù)意義bs表示指定的

8、邊緣線段，如矩形邊界為四段，三角開邊界肯 定為三段。s為第bs段線段弧長(zhǎng)的近似(估計(jì))值，bs與s可以為向量，但是要一一對(duì)應(yīng)，即bs為幾個(gè)值，s也得為幾個(gè)值。輸出變量x,y是每條線段起點(diǎn)和 終點(diǎn)所對(duì)應(yīng)的坐標(biāo)。這個(gè)函數(shù)編制的關(guān)鍵是，函數(shù)內(nèi)邊界上的坐標(biāo)(x(t),y(t)是用參變量t表示的，返回值是求得邊界任意長(zhǎng)度時(shí)的坐標(biāo)(x(t),y(t)值，參量可 以有很多種選法。回到前一篇中，給定的方程的求解區(qū)域是0,1,0,1的一個(gè)正方形，我們將它 命名為mygeom。下面我們來看下 mygeom是怎么編寫的。fun cti on x,y=mygeom(bs,s)nbs=4; %表示邊界的段數(shù)if nar

9、gin=0,x=n bs; %不給定輸入變量時(shí)，輸出表示幾何區(qū)域邊界的線段數(shù)returnendd=000 0 %表示每條線段起始值的參量值t (四條邊界，所以有四列)11 1 1 %表示每條線段的終點(diǎn)值參量值t0 0 0 0 %沿線段方向左邊區(qū)域的標(biāo)識(shí)值，如果是 1,表示選定左邊區(qū)域，如果 %是0,表示不選定左邊區(qū)域1 1 1 1 %沿線段方向右邊區(qū)域的標(biāo)識(shí)值，規(guī)則同上。;bs1= bs(:)'if fin d(bs1<1 | bs1> nbs),error('PDE:squareg:l nv alidBs', 'Non existe nt boun

10、 dary segme nt nu mber.') endif nargin=1,x=d(:,bs1); %給定一個(gè)輸入變量時(shí)，輸出區(qū)域邊界數(shù)據(jù)的矩陣returnendx=zeros(size(s);y=zeros(size(s);m,n =size(bs);if m=1 && n=1,bs=bs* on es(size(s); % expa nd bselseif m=size(s,1) | n=size(s,2),error('PDE:squareg:SizeBs', 'bs must be scalar or of same size as

11、 s.'); endif isempty(s),%第一段邊界ii=fi nd(bs=1);if len gth(ii)x(ii)=i nterp1(d(1,1),d(2,1),0 1,s(ii);% 通過參量來確定邊界上的值y(ii)=i nterp1(d(1,1),d(2,1),1 1,s(ii);%end%第二段邊界ii=fi nd(bs=2);if len gth(ii)x(ii)=i nterp1(d(1,2),d(2,2),1 1,s(ii);y(ii)=in terp1(d(1,2),d(2,2),1 0,s(ii);end%第三段邊界ii=fi nd(bs=3);if l

12、en gth(ii)x(ii)=i nterp1(d(1,3),d(2,3),1 0,s(ii);y(ii)=in terp1(d(1,3),d(2,3),0 0,s(ii);end%第四段邊界ii=fi nd(bs=4);if len gth(ii)x(ii)=i nterp1(d(1,4),d(2,4),0 0,s(ii);y(ii)=in terp1(d(1,4),d(2,4),0 1,s(ii);endend編輯完該函數(shù)后，可以利用pdeplot函數(shù)來測(cè)試設(shè)置的幾何區(qū)域是否滿足要求。在matlab命令窗口輸入：pdegplot('mygeom'), axis equal

13、p,e,t=i nitmesh('mygeom');pdemesh(p,e,t), axis equal結(jié)果如下：00.2DU060.B1mygeom的m文件放在附件里，大家可以下載。本篇寫作時(shí)參考了劉 平的博客，已在前面指出，另參考了偏微分方程的MATLAB軍法及MATLABPDE Toolbox User Guide，有興趣的讀者可以參閱。MATLAB求解PDE問題（3）確定邊界條件已有1508次閱讀2010-10-12 23:46 |個(gè)人分類:生活點(diǎn)滴|系統(tǒng)分類:科研筆記| 關(guān)鍵詞:MATLAB PDE Toolbox pdeboun邊界條件前一篇給出了如何確定幾何區(qū)域，

14、本篇接著給出如何確定邊界條件，在幾何區(qū)域和邊界條件都確定好之后，就可以利用PDE Toolbox對(duì)給定的PDE進(jìn)行計(jì)算了。先回憶下前面的邊界條件：rl O .» / V 9= 0."= i2 far,da ._2L =2-2y- r ,w«*C.Tar上面的四個(gè)邊界條件中，前兩個(gè)是Dirichlet邊界，后兩個(gè)是Neumann邊界 我們先了解下PDE Toolbox中規(guī)定有三種邊界條件，一是Dirichlet邊界條件，一 是廣義Neumann條件，一是混合邊界條件（只用于方程組）。這和傳統(tǒng)的定義 有所不同，按照傳統(tǒng)的觀點(diǎn)，PDE Toolbox中的廣義Neuman

15、n條件應(yīng)該是混合 邊界條件（Robin邊界條件），而傳統(tǒng)的Neumann邊界條件是指邊界處的法向?qū)?數(shù)為零。我們先不考慮方程組，Dirichlet邊界條件和廣義Neumann條件寫成如 下形式：hu = r加（刃其中=心心:血虹門是外汽是邊界法向量與x軸的夾角。規(guī)定邊界條件的m文件是pdebound函數(shù)，該函數(shù)并不提供邊界條件，而是 要求用戶按照規(guī)定的方法給出邊界條件，用戶可以自己給定函數(shù)的名字。我們計(jì)算中給出的邊界條件函數(shù)是 mybound。調(diào)用格式為q,g,h,r=mybound（p,e,u,time）。pdebound函數(shù)中最主要的內(nèi)容是bl矩陣，bl包括了所有的邊界信息。bl的 每一列

16、都是邊界條件矩陣的對(duì)應(yīng)列，每一列都必須滿足下面的規(guī)則：第一行是方程的維數(shù)N，個(gè)方程N(yùn)=1，兩個(gè)方程構(gòu)成的方程組，N=2，； 第二行是Dirichlet邊界條件數(shù) M，Dirichlet邊界條件時(shí)M=N，廣義Neumann 邊界條件時(shí)M=0，如果是方程組，0<MvN之間表示混合邊界條件；第三行到第3+N-1行是表示q的字符串的長(zhǎng)度，這個(gè)長(zhǎng)度按與q有關(guān)的列方 向的次序儲(chǔ)存；2 2第3+N行到第3+N +N-1是表示g的字符串的長(zhǎng)度；2 2第3+N +N行到第3+N +N+MN-1是表示h的字符串的長(zhǎng)度，這個(gè)長(zhǎng)度按與h 有關(guān)的列方向的次序儲(chǔ)存；2 2第3+N +N+MN行到第3+N +N+MN

17、+M-1是表示r的字符串的長(zhǎng)度；接下來的行是包括MATLAB文本表達(dá)式所表示的真實(shí)邊界條件函數(shù)。文本 表達(dá)式是將實(shí)際的表達(dá)式通過 double函數(shù)轉(zhuǎn)化為ASCII碼后的表達(dá)式。文本字 符串的長(zhǎng)度與前面四行規(guī)定的長(zhǎng)度是一致的，兩個(gè)字符串之間沒有分隔符。可以插入含有下列變量的文本表達(dá)式：二維坐標(biāo)x和y;邊界線段參數(shù)s,弧長(zhǎng)的比例。在邊界線段的起始處s=0,向著線段的終點(diǎn)處 增加到s=1;外法向量分量nx, ny。如果要用到切向量，可以用tx和ty表示，其中tx=-ny， ty=nx ；解u （除非已指定輸入?yún)?shù)u）；時(shí)間t （除非已指定輸入?yún)?shù)time）。注意，在只有一個(gè)方程的時(shí)候，即 N=1，如

18、果M=0時(shí)，即是廣義Neumann 邊界條件，此時(shí)不需要表示表示 Dirichlet邊界條件的h，r兩行，因此表示q和 g的字符串從第5行開始（前四行分別為N,M,q,g的指示行）；如果M=1，此時(shí) 表示h，r的字符串從第9行開始（前六行分別是N,M,q,g,h,r的指示行，第七八 行表示q，g，它們均為0，ASCII碼為48）下面看下我們計(jì)算你的例子中邊界mybou nd函數(shù)的內(nèi)容：function q,g,h,r=mybo un d（p,e,u,time） % upper=double('2-2*x-x.A2');% y=1 上邊界條件函數(shù)% down =double(&#

19、39;x.A2'); % y=0 下邊界條件函數(shù)% left =double('y.A2'); % x=0 左邊界條件函數(shù)% right=double('2-2*y-y.A2'); % x=1 右邊界條件函數(shù)% upper=upper' %轉(zhuǎn)化為列向量% dow n =dow n'% left =left'% right=right'%bl=1 1 1 1%0 0 1 1%界條件1 1 1 1%10 10 1 1% g48 48 1 1% h50 50 4 4% r45 45 48 48%上面50 50 48 48%數(shù)42

20、 42 49 49%120121120121 行開始，45 45 46 46即7,8兩行是120 121 94 94的 ASCII46 46 50 50N表示方程維數(shù)的指示行M 表示Dirichlet邊界條件的指示行，M=0為廣義Neumann邊q表示q的字符串的長(zhǎng)度，根據(jù)實(shí)際長(zhǎng)度確定，表示g的字符串的長(zhǎng)度，如1表示字符串的長(zhǎng)度為1表示h的字符串的長(zhǎng)度，這是針對(duì) M=1的列，本例子是三四列 表示r的字符串的長(zhǎng)度，規(guī)則同上。本行以下是邊界條件函數(shù)的字符串表達(dá)式，字符串的長(zhǎng)度與規(guī)定的一致。注意，這些都是將表達(dá)式轉(zhuǎn)化為 ASCII碼后的例如x的ASCII碼是120, ASCII碼的0表示null （

21、空）%需要注意的還有，M=0的列，字符表達(dá)式從g的下一%即從第5行開始。M=1的列，從r的下一行開始，% 表示q，g的字符串的長(zhǎng)度，此時(shí)認(rèn)為q，g均為0，0%碼是48,因此7,8兩行均為48。在往下，按照h，r的長(zhǎng)度排列。94 94 0 0齊，一般以050 50 0 0%由此造成的行數(shù)不匹配，可以由任意 ASCII碼補(bǔ)% （表示null,空）表示，不容易引起歧義。； if any （size（u）q,g,h,r=pdeexpd（p,e,u,time,bl）;elseq,g,h,r=pdeexpd（p,e,time,bl）;end規(guī)定邊界條件的矩陣bl較為復(fù)雜，詳細(xì)的解釋可以使大家省不少時(shí)間，

22、但是 也容易把大家整暈，我也是在這方面花費(fèi)了不少功夫。下面具體分析下bl中的兩列。第一列對(duì)應(yīng)著上邊界條件（Neuma nn邊界條件）：MATLAB求解PDE問題（4）總結(jié)與邊界條件的標(biāo)準(zhǔn)形式比較，有，q=0, g=2-2*x-x.A2，于是第一列可以寫為1 0 1 10 '0''2-2*x-x.A3''其中1表示N=1，即一個(gè)方程，0表示Neumann邊界條件；1表示q的長(zhǎng)度；10 表示g的長(zhǎng)度，0表示q=0，長(zhǎng)度為1； '2-2*x-x.A3'表示g=2-2*x-x.A2，轉(zhuǎn)化為ASCII碼后長(zhǎng)度為10。因此第一列最終寫為1 0 1 1

23、0 48 50 45 50 42 120 45 120 46 94 50 '第三列對(duì)應(yīng)著下邊界條件，是 Dirichlet邊界條件：與標(biāo)準(zhǔn)形式比較，有，h=1，r=x.A2，于是第三列可以寫為1 1 1 1 1 4 '0' '0' '1' 'x.A2''按順序依次為，1表示N=1，個(gè)方程；1表示Dirichlet邊界條件；1表示q的 長(zhǎng)度，此時(shí)必為1； 1表示g的長(zhǎng)度，此時(shí)也必為1； 1表示h的長(zhǎng)度是1； 4表 示r的長(zhǎng)度為4； '0'表示q的值為0 （0的字符串長(zhǎng)度是1）； '0'

24、表示g的值是0；'1'表示h=1（ 1的ASCII碼是49）； 'x.A2'表示r=x.A2，轉(zhuǎn)化為ASCII碼后長(zhǎng)度是4。于是第三列可以寫為1 1 1 1 1 4 48 48 49 120 46 94 50'此時(shí)第三列的長(zhǎng)度小于第一列的長(zhǎng)度，缺幾個(gè)長(zhǎng)度就可以添加幾個(gè) ASCII碼的0, 表示空，這樣第三列就是 mybound中的第三列了，如下1 1 1 1 1 4 48 48 49 120 46 94 50 0 0'最后附上mybound的m文件，有需要的可以下載。mybo und 文件MATLAB求解PDE問題（4）總結(jié)已有1792次閱讀20

25、10-10-13 11:16 |個(gè)人分類:生活點(diǎn)滴|系統(tǒng)分類:科研筆記 | 關(guān)鍵詞:MATLAB PDE Toolbo使用說明 assempde通過前面三篇的介紹，大家已經(jīng)能夠掌握PDE Toolbox求解橢圓型方程的用 法了。本篇在前面三篇的基礎(chǔ)上，大致介紹PDE Toolbox中一些其它函數(shù)的用法， 供大家理解。我們先總結(jié)下利用PDE Toolbox求解PDE(s)的基本步驟：第一要把方程化為 標(biāo)準(zhǔn)形式，不論是橢圓方程，還是其他類型的方程，均要化為標(biāo)準(zhǔn)形式后才能利 用工具箱求解；第二要寫出表示求解區(qū)域和幾何區(qū)域文件以及表示邊界條件的邊 界條件文件；完成這兩步后，后面的步驟就容易了。依次為網(wǎng)

26、格初始化，加密網(wǎng) 格(或者使用自適應(yīng)網(wǎng)格)；根據(jù)方程的類型使用相應(yīng)的計(jì)算函數(shù)；對(duì)結(jié)果進(jìn)行 分析計(jì)算以及可視化處理。PDE Toolbox能夠求解的方程類型前面也已說過，有橢圓型方程，拋物型方 程，雙曲型方程，特征值方程，非線性方程五類。其中拋物型和雙曲型涉及到時(shí) 間，應(yīng)按照要求給出初始值。特征值問題是一個(gè)齊次問題，即邊界條件中g(shù)=0,r=0。非齊次部分將會(huì)被自動(dòng)刪去，這在理解工具箱自帶的squareb2邊界條件時(shí)尤為重要。下面以一段程序來盡量多介紹幾個(gè)函數(shù)的用法，關(guān)于這些函數(shù)的用法大家可以從網(wǎng)絡(luò)、工具箱的用戶說明手冊(cè)、各種介紹工具箱的數(shù)據(jù)中得到詳細(xì)的說明。g='yourgeom'b='yourbo un d'以上兩行確定幾何區(qū)域和邊界條件p,e,t=i ni tmesh('yourgeom');% 網(wǎng)格初始化p,e,t= refinemesh('yourgeom',p,e,t);% 網(wǎng)格加密p,e,t= refin emesh('yourgeom',p,e,t);以上三行初始化網(wǎng)格并加密c='yourcontent c' %系數(shù)a='yourcontent a'f='yourcontent f;以上三行確定方程參數(shù)，注意特征值問題是齊次問題，不需要