應(yīng)用Matlab的PDE Toolbox求解偏微分方程綜合實(shí)驗(yàn)

1、2011屆信計(jì)專業(yè)學(xué)生綜合實(shí)驗(yàn)題目（要求按照所附開題報(bào)告表的格式填寫提交開題報(bào)告。一個(gè)小組選做一題，小組全體成員共同完成，每個(gè)小組只提交一份實(shí)驗(yàn)報(bào)告，按照出力多少排名。提交時(shí)間在本學(xué)期18周以前。）3 應(yīng)用matlab的pde toolbox求解偏微分方程熟悉matlab的pde工具箱的功能，并用其求解具有工程背景的偏微分方程，要求分別對(duì)三種類型方程：拋物型、橢圓形和雙曲型。闡述清楚如下方法： 1、這里我們先脫離問題所含有的工程背景，分別舉三個(gè)例子，大致的說明一下如何使用matlab的pde工具箱來求解三種類型的偏微分方程。(1) 橢圓形方程：考慮一塊圓形金屬片，中心挖去一正方形，外邊界滿足n

2、eumann條件，內(nèi)邊界滿足dirichlet條件：考慮到入射波以方向，所以上式可以寫成這樣得到求解這個(gè)入射波的定解問題：這里取波長為0.1。 現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形： 圖1 初始網(wǎng)格、加密以及網(wǎng)格剖分?jǐn)?shù)據(jù) 圖2 解的三維圖形 圖3 解的二維動(dòng)畫圖 附錄1：二維動(dòng)畫的matlab程序echo onclc%程序段一：解helmholtz方程-div(grad(u)-k2u=0%并研究正方形上的反射波，波源來自右邊clc%入射波波數(shù)為60k=60;g=scatterg;%scatterg描述幾何區(qū)域的文件名，此區(qū)域?yàn)閳@內(nèi)有一方洞b=scatterb;%scatterb

3、是描述邊界條件的文件名，內(nèi)邊界滿足 %dirichlet條件，外邊界滿足neumann條件%選擇方程系數(shù)c,a,fc=1;a=-k2;f=0;%程序段二：初始化網(wǎng)格和加密網(wǎng)格p,e,t=initmesh(g);p,e,t=refinemesh(g,p,e,t);p,e,t=refinemesh(g,p,e,t);%會(huì)出網(wǎng)格圖pdemesh(p,e,t);axis equalclc%在復(fù)平面上求解u=assempde(b,p,e,t,c,a,f);%取復(fù)數(shù)解的實(shí)部h=newplot;set(get(h,parent),renderer,zbuffer)pdeplot(p,e,t,xydata,r

4、eal(u),zdata,real(u),. mesh,off);colormap(cool)clc%制作反射波的動(dòng)畫程序m=10;%幀數(shù)h=newplot;hf=get(h,parent);set(hf,renderer,zbuffer)axis tight,set(gca,dataaspectratio,1 1 1);axis offm=moviein(m,hf);maxu=max(abs(u);for j=1:m,. uu=real(exp(-j*2*pi/m*sqrt(-1)*u);. fprintf(%d,j);. pdeplot(p,e,t,xydata,uu,colorbar,o

5、ff,. mesh,off),. caxis(-maxu maxu);. axis tight,set(gca,dataaspectratio,1 1 1);. axis off,. m(:,j)=getframe(hf);. if j=m,. fprintf(donen);. end,.end%顯示動(dòng)畫movie(hf,m,50);echo off(2) 拋物型方程：考慮一個(gè)圓柱形放射性桿，其左端供熱，右端保持常溫，側(cè)面與環(huán)境有熱交換。由于放性作用，熱量均勻地產(chǎn)生。初始溫度為。于是可以用如下方程描述： 其中為密度，為桿的熱容量，為導(dǎo)熱系數(shù)，為放射性熱源密度。這一金屬桿的密度 取為熱容量為導(dǎo)熱

6、系數(shù)為熱源密度為右端恒溫為側(cè)面環(huán)境溫度為熱交換系數(shù)為左端的熱流為 邊界條件(如右圖4)： 在桿的左端(): 在桿的右端(): 在桿的側(cè)面()： 在桿的軸心(): 初始溫度： 現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形：圖5 解的2維動(dòng)畫圖形 圖6 解的3維動(dòng)畫圖形 圖7 解的動(dòng)畫圖形比較圖(3) 雙曲型方程：考慮如下二維波動(dòng)方程的定界問題，并最終獲得其解得圖形：現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形：圖8 解的3維動(dòng)畫圖形附錄2：雙曲線三維動(dòng)畫的matlab程序p.e.t=initmesh(squareg);%初始化網(wǎng)格x=p(1,:),y=p(2,:);u0=atan

7、(cos(pi/2*x);%ut0=3*sin(pi*x).*exp(sin(pi/2*y);%n=31;tlist=linspace(0,5,n);%uu=hyperbolic(u0,ut0,tlist,squareb3,p,e,t,1,0,0,1);%delta=-1:0.1:1;uxy,tn,a2,a3=tri2grid(p,t,uu(:,1),detla,delta);%gp=tn;a2;a3;umax=max(max(uu);umin=min(min(uu);newplotm=moviein(n);for i=1:n, pdeplot(p,e,t,xydata,uu(:,i),zda

8、ta,uu(:,i),mesh,off,xygrid,on,gridparam,gp,colorbar,off,zstyle,continuous);axis(-1 1 -1 1 umin umax);caxis(umin umax);m(:,i)=getframe;endmovie(m,10);2、向問題中添加一定的工程背景，以下我們?cè)敱M的說明如何使用matlab的pde工具箱來求解三種類型的偏微分方程。(1)橢圓形方程： 在磁鐵、電動(dòng)機(jī)、變壓器這類問題中，考慮如下定常的maxwell方程：其中b是磁感應(yīng)強(qiáng)度，h是磁場(chǎng)強(qiáng)度，j是電流密度，是材料的磁導(dǎo)率。由于故存在一個(gè)靜磁矢勢(shì)a，使得 及 平

9、面問題中假設(shè)電流平行于z軸，故a僅有z分量 從而，上述方程可以簡(jiǎn)化為橢圓型pde： 其中 對(duì)于二維情況，可以計(jì)算磁感應(yīng)強(qiáng)度為磁場(chǎng)強(qiáng)度為 對(duì)于不同性質(zhì)的材料組成的區(qū)域之間的交界面上，是連續(xù)的，也就是連續(xù)。 在強(qiáng)磁場(chǎng)性材料中，依賴于場(chǎng),從而導(dǎo)致非線性解。 如果是dirichlet邊界條件，則要求在邊界上給出磁勢(shì)a的值。如果是neumannb邊界條件，則要求在邊界上給出的值，也就是在邊界上給出磁場(chǎng)強(qiáng)度h的切向分量。磁勢(shì)a、磁場(chǎng)強(qiáng)度h以及磁感應(yīng)強(qiáng)度b的可視化均可實(shí)現(xiàn)，還可以作b和h的向量場(chǎng)圖形。 考慮由雙極電動(dòng)錠子線圈產(chǎn)生的靜磁場(chǎng)問題。假設(shè)馬達(dá)很長，使得端點(diǎn)的影響可以忽略不計(jì)，從而簡(jiǎn)化成二維模型處理。

10、 區(qū)域由四部分組成： 兩個(gè)鐵磁體：定子和轉(zhuǎn)子 定子和轉(zhuǎn)子之間的間隙 帶有直流的轉(zhuǎn)子線圈在空氣和線圈中磁導(dǎo)率為1，而定子和轉(zhuǎn)子的磁導(dǎo)率為其中這些是由變壓器鋼材決定。 線圈的電流密度j為1，其他均為0。這一問題比較復(fù)雜，它的區(qū)域包括5個(gè)圓和兩個(gè)矩形(如圖9所示)，它的方程在不同區(qū)域也不完全相同，而且包括非線性方程。以下我們?cè)敿?xì)說明利用gui求解的過程，依次回答以下問題：圖9(1) 如何啟動(dòng)pde toolbox工具箱： 啟動(dòng)matlab7.0，在其命令窗口(command window)中鍵入命令pdetool，這時(shí)系統(tǒng)將立即產(chǎn)生偏微分方程工具箱(pde toolbox)的圖形用戶界面(graph

11、ical user interface,簡(jiǎn)記為gui)，也就實(shí)現(xiàn)了pde toolbox工具箱的啟動(dòng)。 此外，由于此問題所涉及的是靜磁場(chǎng)的問題，因此，我們需要將系統(tǒng)默認(rèn)的應(yīng)用模型模式選定為magnetostatics.(2)怎樣分步完成平面幾何造型： 由于該問題所涉及的區(qū)域方程包括非線性方程，因此我們得選擇使用非線性求解，即單擊solve菜單中paraneters.選項(xiàng),打開solve parameters對(duì)話框，選定use nonlinear solver. 在pde toolbox面板中依次拖出圓雙擊所繪制的圓，將圓心設(shè)置為(0,0),半徑分別為1,0.8,0.6,0.5,0.4；依次拖出

12、矩形雙擊所繪制的矩形，分別設(shè)置left為-0.2、-0.1，bottom分別為0.2、0.1,width分別為0.2、0.2，height分別為0.9、0.9.這樣便可以得到圖9所示的平面圖. 此外上述圖形繪制過程還可以通過向matlab中鍵入如下命令實(shí)現(xiàn)： 附錄3pdecirc(0,0,1,c1)pdecirc(0,0,0.8,c2)pdecirc(0,0,0.6,c3)pdecirc(0,0,0.5,c4)pdecirc(0,0,0.4,c5)pderect(-0.2 0.2 0.2 0.9,r1)pderect(-0.1 0.1 0.2 0.9,r2)pderect(0 1 0 1,sq

13、1) 單擊boundary菜單下的boundary mode選項(xiàng)，進(jìn)入邊界模式.選擇boundary菜單中的show subdomain labels命令，顯示區(qū)域編號(hào)(如圖10所示)，再單擊所要?jiǎng)h去的線段或弧段，執(zhí)行boundary菜單中remove subdomain border命令，則可以得到圖11. 圖10圖11 為了便于圖形的顯示，為此我們重新調(diào)整坐標(biāo)軸的范圍，單擊options菜單中axes limits. 選項(xiàng)，打開axes limits對(duì)話框，取消auto選項(xiàng)，取為默認(rèn)值。再單擊options菜單中axes equal選項(xiàng)，這樣所得到的兩坐標(biāo)軸的尺寸比例一致。(3)如何選取邊

14、界： 逐段雙擊邊界設(shè)置邊界條件：x軸邊界選齊次neumann條件，鍵入其他邊界選齊次dirichlet條件，系統(tǒng)已經(jīng)默認(rèn)這種邊界條件.(4) 如何進(jìn)入pde模式：用gui設(shè)置定解問題包括三個(gè)模式(mode):(1) draw模式：使用csg(幾何結(jié)構(gòu)實(shí)體模型)對(duì)話框畫幾何區(qū)域，包括矩形、圓、橢圓和多邊形，也可以將它們組合使用。(2) boundary模式：在各個(gè)邊界段上給出邊界條件。(3) pde模式：確定方程的類型、c,a,f和d.也能夠在不同的子區(qū)域上 設(shè)置不同的系數(shù)(反映材料的性質(zhì)).用gui設(shè)置定解問題包括三個(gè)模式(mode): (1)mesh模式：生成網(wǎng)格、自動(dòng)控制參數(shù)。 (2)so

15、lve模式：對(duì)于橢圓型方程還能求非線性和自適應(yīng)解。對(duì)于拋物型和雙曲型方程，設(shè)置初始邊值條件后能求出給定t時(shí)刻的解。求解后可以加密網(wǎng)格再求解。為此，選擇pde菜單中pde mode命令，進(jìn)入pde模式(如圖12)，此時(shí)需要確定方程的類型、c,a,f和d.也能夠在不同的子區(qū)域上設(shè)置不同的系數(shù)(反映材料的性質(zhì))：區(qū)域4鍵入mu=1,j=1;區(qū)域1、2均鍵入這里系數(shù)含有非線性；區(qū)域3鍵入mu=1,j=0.圖12 (5)如何初始化網(wǎng)格： 選擇mesh菜單中initialize mesh命令，即可對(duì)平面圖形進(jìn)行網(wǎng)格剖分. (6)如何進(jìn)行網(wǎng)格加密： 選擇mesh菜單中refine mesh命令，即可實(shí)現(xiàn)對(duì)網(wǎng)

16、格進(jìn)行加密. (7)怎樣進(jìn)行求解，如何顯示圖形： 選擇solve菜單中solve pde命令，解偏微分方程并顯示圖形解(如圖13所示),圖13(8)顯示三維圖形解： 單擊plot菜單中parameters.選項(xiàng)，打開plot selection對(duì)話框，選中color,contour ,height(3-d plot)和arrows四項(xiàng)，設(shè)置contour plot levels:10,然后單擊plot按鈕，顯示三維圖形解.如圖14所示圖14 解的3維圖形(2)拋物型方程：考慮拋物線偏微分方程：它描述平面的熱傳輸現(xiàn)象，以及由軸對(duì)稱三維問題經(jīng)過降維后的熱傳輸問題，其中t為溫度，其他參數(shù)為： 在pd

17、e toolbox中我們可以做溫度、溫度梯度、熱流的可視化圖形。 考慮一個(gè)方形區(qū)域，導(dǎo)熱系數(shù)為10，密度為2。在方形區(qū)域內(nèi)有一個(gè)圓形區(qū)域的熱源為4，導(dǎo)熱系數(shù)為2，密度為1。兩個(gè)區(qū)域的比熱都為0.1。以下我們?cè)敿?xì)說明利用gui求解的過程，依次回答以下問題：(1)如何啟動(dòng)pde toolbox工具箱：?jiǎn)?dòng)matlab7.0，在其命令窗口(command window)中鍵入命令pdetool，這時(shí)系統(tǒng)將立即產(chǎn)生偏微分方程工具箱(pde toolbox)的圖形用戶界面(graphical user interface,簡(jiǎn)記為gui)，也就實(shí)現(xiàn)了pde toolbox工具箱的啟動(dòng)。此外，由于此問題所涉

18、及的是平面熱傳輸?shù)膯栴}，因此，我們需要將系統(tǒng)默認(rèn)的應(yīng)用模型模式選定為heat transfer. (2)怎樣分步完成平面幾何造型： 為了便于平面圖形的觀察，為此，我們對(duì)坐標(biāo)的范圍顯示進(jìn)行設(shè)置，即單擊options菜單中的grid spacing對(duì)話框中選auto;打開axes limits對(duì)話框，設(shè)置x軸和y軸的范圍為0,3,單擊apply即可完成設(shè)置。 此外,為了使兩坐標(biāo)軸的尺寸比例一致,我們重新調(diào)整坐標(biāo)軸的范圍，即單擊options菜單中axes equal選項(xiàng)。 在pde toolbox面板中依次拖出圓方形分別雙擊所繪制的圖形，對(duì)于圓:設(shè)置x-center為1.5，y-center為1.

19、5，radius為1;對(duì)于方形;設(shè)置left為0，bottom分別為0,width分別為3，height分別為3. 這樣便可以得到圖15所示的平面圖.圖15 因?yàn)樗型膺吔绲臏囟缺３譃?，可知邊界滿足默認(rèn)的齊次dirichlet邊界條件。所以,只需單擊boundary菜單下的boundary mode選項(xiàng)，進(jìn)入邊界模式.采取系統(tǒng)默認(rèn)的邊界方式。 在pde mode模式下，分別雙擊兩個(gè)區(qū)域設(shè)置拋物型pde參數(shù):依次將c,a,f,d,設(shè)為2,0,4,0.1;10,0,0,0.2. (3)如何選取邊界： 因?yàn)樗型膺吔绲臏囟缺３譃?，可知邊界滿足默認(rèn)的齊次dirichlet邊界條件。只需單擊boun

20、dary菜單下的boundary mode選項(xiàng)，進(jìn)入邊界模式.這里我們采取系統(tǒng)默認(rèn)的邊界方式。 (4)如何進(jìn)入pde模式：用gui設(shè)置定解問題包括三個(gè)模式(mode):(4) draw模式：使用csg(幾何結(jié)構(gòu)實(shí)體模型)對(duì)話框畫幾何區(qū)域，包括矩形、圓、橢圓和多邊形，也可以將它們組合使用。(5) boundary模式：在各個(gè)邊界段上給出邊界條件。(6) pde模式：確定方程的類型、c,a,f和d.也能夠在不同的子區(qū)域上 設(shè)置不同的系數(shù)(反映材料的性質(zhì)).用gui設(shè)置定解問題包括三個(gè)模式(mode): (1)mesh模式：生成網(wǎng)格、自動(dòng)控制參數(shù)。 (2)solve模式：對(duì)于橢圓型方程還能求非線性和

21、自適應(yīng)解。對(duì)于拋物型和雙曲型方程，設(shè)置初始邊值條件后能求出給定t時(shí)刻的解。求解后可以加密網(wǎng)格再求解。為此，選擇pde菜單中pde mode命令，進(jìn)入pde模式，此時(shí)需要確定方程的類型、c,a,f和d. (5)如何初始化網(wǎng)格： 選擇mesh菜單中initialize mesh命令，即可對(duì)平面圖形進(jìn)行網(wǎng)格剖分. (6)如何進(jìn)行網(wǎng)格加密： 選擇mesh菜單中refine mesh命令，即可實(shí)現(xiàn)對(duì)網(wǎng)格進(jìn)行加密. (7)怎樣進(jìn)行求解，如何顯示圖形： 選擇solve菜單中solve pde命令，解偏微分方程并顯示圖形解(如圖16所示), 圖16(8)顯示三維圖形解： 單擊plot菜單中parameters

22、.選項(xiàng)，打開plot selection對(duì)話框，選中color,height(3-d plot)、arrows和animation四項(xiàng)，設(shè)置contour plot levels:10,然后單擊plot按鈕，顯示三維圖形動(dòng)畫解.如圖17所示 圖17 解的三維圖形動(dòng)畫解(3)雙曲型方程： 此例可參見第一種情形第三種情況,以下8個(gè)單問均可參考第二種情形兩種情況,這里作省略. 云南財(cái)經(jīng)大學(xué)實(shí)驗(yàn)報(bào)告系 （院）： 統(tǒng)數(shù)學(xué)院 專 業(yè)： 班 級(jí)： 姓 名： 課程名稱： 復(fù)合實(shí)驗(yàn) 實(shí)驗(yàn)時(shí)間： 2010-12-12 指導(dǎo)教師： 云南財(cái)經(jīng)大學(xué)教務(wù)處制填表說明1 實(shí)驗(yàn)名稱 要用最簡(jiǎn)練的語言反映實(shí)驗(yàn)的內(nèi)容。2 實(shí)驗(yàn)?zāi)?/p>

23、的 目的要明確，要抓住重點(diǎn)，可以從理論和實(shí)踐兩個(gè)方面考慮。在理論上，驗(yàn)證定理、公式、算法，并使實(shí)驗(yàn)者獲得深刻和系統(tǒng)的理解，在實(shí)踐上，掌握使用實(shí)驗(yàn)設(shè)備的技能技巧和程序調(diào)試的方法。一般需要說明是驗(yàn)證型實(shí)驗(yàn)還是設(shè)計(jì)型實(shí)驗(yàn)，是創(chuàng)新型實(shí)驗(yàn)還是綜合型實(shí)驗(yàn)。3 實(shí)驗(yàn)環(huán)境 實(shí)驗(yàn)用的軟硬件環(huán)境（配置）。4 實(shí)驗(yàn)內(nèi)容（算法、程序、步驟和方法）這是實(shí)驗(yàn)報(bào)告極其重要的內(nèi)容。這部分要寫明依據(jù)何種原理、定律算法、或操作方法進(jìn)行實(shí)驗(yàn)，要寫明經(jīng)過哪幾個(gè)步驟。還應(yīng)該畫出流程圖（實(shí)驗(yàn)裝置的結(jié)構(gòu)示意圖），再配以相應(yīng)的文字說明，這樣既可以節(jié)省許多文字說明，又能使實(shí)驗(yàn)報(bào)告簡(jiǎn)明扼要，清楚明白。5 結(jié)論（結(jié)果） 即根據(jù)實(shí)驗(yàn)過程中所見到現(xiàn)象

24、和測(cè)得的數(shù)據(jù)，作出結(jié)論。6 小結(jié) 對(duì)本次實(shí)驗(yàn)的思考和建議。7 備注或說明 可填寫實(shí)驗(yàn)成功或失敗的原因，實(shí)驗(yàn)后的心得體會(huì)等。8 指導(dǎo)教師評(píng)分 指導(dǎo)教師根據(jù)本次實(shí)驗(yàn)的預(yù)習(xí)、表現(xiàn)、操作的實(shí)驗(yàn)報(bào)告的撰寫客觀進(jìn)行評(píng)分、簽名，并記入成績(jī)。實(shí)驗(yàn)名稱應(yīng)用matlab的pde toolbox求解偏微分方程實(shí)驗(yàn)?zāi)康氖煜atlab的pde工具箱的功能，并用其求解具有工程背景的偏微分方程，要求分別對(duì)三種類型方程：拋物型、橢圓形和雙曲型。闡述清楚如下方法：(1)如何啟動(dòng)pde toolbox工具箱：(2)怎樣分步完成平面幾何造型：(3)如何選取邊界：(4)如何進(jìn)入pde模式：(5)如何初始化網(wǎng)格：(6)如何進(jìn)行網(wǎng)格加

25、密：(7)怎樣進(jìn)行求解，如何顯示圖形：(8)顯示三維圖形解：實(shí)驗(yàn)內(nèi)容（算法、程序、步驟和方法）1、這里我們先脫離問題所含有的工程背景，分別舉三個(gè)例子，大致的說明一下如何使用matlab的pde工具箱來求解三種類型的偏微分方程。(1)橢圓形方程：考慮一塊圓形金屬片，中心挖去一正方形，外邊界滿足neumann條件，內(nèi)邊界滿足dirichlet條件：考慮到入射波以方向，所以上式可以寫成這樣得到求解這個(gè)入射波的定解問題：這里取波長為0.1。 現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形： 圖1 初始網(wǎng)格、加密以及網(wǎng)格剖分?jǐn)?shù)據(jù) 圖2 解的三維圖形 圖3 解的二維動(dòng)畫圖 附錄1：二維動(dòng)畫的matl

26、ab程序echo onclc%程序段一：解helmholtz方程-div(grad(u)-k2u=0%并研究正方形上的反射波，波源來自右邊clc%入射波波數(shù)為60k=60;g=scatterg;%scatterg描述幾何區(qū)域的文件名，此區(qū)域?yàn)閳@內(nèi)有一方洞b=scatterb;%scatterb是描述邊界條件的文件名，內(nèi)邊界滿足 %dirichlet條件，外邊界滿足neumann條件%選擇方程系數(shù)c,a,fc=1;a=-k2;f=0;%程序段二：初始化網(wǎng)格和加密網(wǎng)格p,e,t=initmesh(g);p,e,t=refinemesh(g,p,e,t);p,e,t=refinemesh(g,p,e

27、,t);%會(huì)出網(wǎng)格圖pdemesh(p,e,t);axis equalclc%在復(fù)平面上求解u=assempde(b,p,e,t,c,a,f);%取復(fù)數(shù)解的實(shí)部h=newplot;set(get(h,parent),renderer,zbuffer)pdeplot(p,e,t,xydata,real(u),zdata,real(u),. mesh,off);colormap(cool)clc%制作反射波的動(dòng)畫程序m=10;%幀數(shù)h=newplot;hf=get(h,parent);set(hf,renderer,zbuffer)axis tight,set(gca,dataaspectrati

28、o,1 1 1);axis offm=moviein(m,hf);maxu=max(abs(u);for j=1:m,. uu=real(exp(-j*2*pi/m*sqrt(-1)*u);. fprintf(%d,j);. pdeplot(p,e,t,xydata,uu,colorbar,off,. mesh,off),. caxis(-maxu maxu);. axis tight,set(gca,dataaspectratio,1 1 1);. axis off,. m(:,j)=getframe(hf);. if j=m,. fprintf(donen);. end,.end%顯示動(dòng)畫

29、movie(hf,m,50);echo off(2)拋物型方程：考慮一個(gè)圓柱形放射性桿，其左端供熱，右端保持常溫，側(cè)面與環(huán)境有熱交換。由于放性作用，熱量均勻地產(chǎn)生。初始溫度為。于是可以用如下方程描述： 其中為密度，為桿的熱容量，為導(dǎo)熱系數(shù)，為放射性熱源密度。這一金屬桿的密度 取為熱容量為導(dǎo)熱系數(shù)為熱源密度為右端恒溫為側(cè)面環(huán)境溫度為熱交換系數(shù)為左端的熱流為 邊界條件(如右圖4)： 在桿的左端(): 在桿的右端(): 在桿的側(cè)面()： 在桿的軸心(): 初始溫度： 現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形：圖5 解的2維動(dòng)畫圖形 圖6 解的3維動(dòng)畫圖形 圖7 解的動(dòng)畫圖形比較圖(3)雙

30、曲型方程：考慮如下二維波動(dòng)方程的定界問題，并最終獲得其解得圖形：現(xiàn)在用gui來求解上述這個(gè)問題，并最終獲得其解得圖形：圖8 解的3維動(dòng)畫圖形附錄2：雙曲線三維動(dòng)畫的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(sin(pi/2*y);%n=31;tlist=linspace(0,5,n);%uu=hyperbolic(u0,ut0,tlist,squareb3,p,e,t,1,0,0,1);%delta=-1:0.1:1;uxy,tn,a2,a3

31、=tri2grid(p,t,uu(:,1),detla,delta);%gp=tn;a2;a3;umax=max(max(uu);umin=min(min(uu);newplotm=moviein(n);for i=1:n, pdeplot(p,e,t,xydata,uu(:,i),zdata,uu(:,i),mesh,off,xygrid,on,gridparam,gp,colorbar,off,zstyle,continuous);axis(-1 1 -1 1 umin umax);caxis(umin umax);m(:,i)=getframe;endmovie(m,10);2、向問題中

32、添加一定的工程背景，以下我們?cè)敱M的說明如何使用matlab的pde工具箱來求解三種類型的偏微分方程。(1)橢圓形方程： 在磁鐵、電動(dòng)機(jī)、變壓器這類問題中，考慮如下定常的maxwell方程：其中b是磁感應(yīng)強(qiáng)度，h是磁場(chǎng)強(qiáng)度，j是電流密度，是材料的磁導(dǎo)率。由于故存在一個(gè)靜磁矢勢(shì)a，使得 及 平面問題中假設(shè)電流平行于z軸，故a僅有z分量 從而，上述方程可以簡(jiǎn)化為橢圓型pde： 其中 對(duì)于二維情況，可以計(jì)算磁感應(yīng)強(qiáng)度為磁場(chǎng)強(qiáng)度為 對(duì)于不同性質(zhì)的材料組成的區(qū)域之間的交界面上，是連續(xù)的，也就是連續(xù)。 在強(qiáng)磁場(chǎng)性材料中，依賴于場(chǎng),從而導(dǎo)致非線性解。 如果是dirichlet邊界條件，則要求在邊界上給出磁勢(shì)a

33、的值。如果是neumannb邊界條件，則要求在邊界上給出的值，也就是在邊界上給出磁場(chǎng)強(qiáng)度h的切向分量。磁勢(shì)a、磁場(chǎng)強(qiáng)度h以及磁感應(yīng)強(qiáng)度b的可視化均可實(shí)現(xiàn)，還可以作b和h的向量場(chǎng)圖形。 考慮由雙極電動(dòng)錠子線圈產(chǎn)生的靜磁場(chǎng)問題。假設(shè)馬達(dá)很長，使得端點(diǎn)的影響可以忽略不計(jì)，從而簡(jiǎn)化成二維模型處理。 區(qū)域由四部分組成： 兩個(gè)鐵磁體：定子和轉(zhuǎn)子 定子和轉(zhuǎn)子之間的間隙 帶有直流的轉(zhuǎn)子線圈在空氣和線圈中磁導(dǎo)率為1，而定子和轉(zhuǎn)子的磁導(dǎo)率為其中這些是由變壓器鋼材決定。 線圈的電流密度j為1，其他均為0。這一問題比較復(fù)雜，它的區(qū)域包括5個(gè)圓和兩個(gè)矩形(如圖9所示)，它的方程在不同區(qū)域也不完全相同，而且包括非線性方程

34、。以下我們?cè)敿?xì)說明利用gui求解的過程，依次回答以下問題：圖9(2) 如何啟動(dòng)pde toolbox工具箱： 啟動(dòng)matlab7.0，在其命令窗口(command window)中鍵入命令pdetool，這時(shí)系統(tǒng)將立即產(chǎn)生偏微分方程工具箱(pde toolbox)的圖形用戶界面(graphical user interface,簡(jiǎn)記為gui)，也就實(shí)現(xiàn)了pde toolbox工具箱的啟動(dòng)。 此外，由于此問題所涉及的是靜磁場(chǎng)的問題，因此，我們需要將系統(tǒng)默認(rèn)的應(yīng)用模型模式選定為magnetostatics.(2)怎樣分步完成平面幾何造型： 由于該問題所涉及的區(qū)域方程包括非線性方程，因此我們得選擇使

35、用非線性求解，即單擊solve菜單中paraneters.選項(xiàng),打開solve parameters對(duì)話框，選定use nonlinear solver. 在pde toolbox面板中依次拖出圓雙擊所繪制的圓，將圓心設(shè)置為(0,0),半徑分別為1,0.8,0.6,0.5,0.4；依次拖出矩形雙擊所繪制的矩形，分別設(shè)置left為-0.2、-0.1，bottom分別為0.2、0.1,width分別為0.2、0.2，height分別為0.9、0.9.這樣便可以得到圖9所示的平面圖. 此外上述圖形繪制過程還可以通過向matlab中鍵入如下命令實(shí)現(xiàn)： 附錄3pdecirc(0,0,1,c1)pdeci

36、rc(0,0,0.8,c2)pdecirc(0,0,0.6,c3)pdecirc(0,0,0.5,c4)pdecirc(0,0,0.4,c5)pderect(-0.2 0.2 0.2 0.9,r1)pderect(-0.1 0.1 0.2 0.9,r2)pderect(0 1 0 1,sq1) 單擊boundary菜單下的boundary mode選項(xiàng)，進(jìn)入邊界模式.選擇boundary菜單中的show subdomain labels命令，顯示區(qū)域編號(hào)(如圖10所示)，再單擊所要?jiǎng)h去的線段或弧段，執(zhí)行boundary菜單中remove subdomain border命令，則可以得到圖11.

37、 圖10圖11 為了便于圖形的顯示，為此我們重新調(diào)整坐標(biāo)軸的范圍，單擊options菜單中axes limits. 選項(xiàng)，打開axes limits對(duì)話框，取消auto選項(xiàng)，取為默認(rèn)值。再單擊options菜單中axes equal選項(xiàng)，這樣所得到的兩坐標(biāo)軸的尺寸比例一致。(3)如何選取邊界： 逐段雙擊邊界設(shè)置邊界條件：x軸邊界選齊次neumann條件，鍵入其他邊界選齊次dirichlet條件，系統(tǒng)已經(jīng)默認(rèn)這種邊界條件.(4)如何進(jìn)入pde模式：用gui設(shè)置定解問題包括三個(gè)模式(mode):(1)draw模式：使用csg(幾何結(jié)構(gòu)實(shí)體模型)對(duì)話框畫幾何區(qū)域，包括矩形、圓、橢圓和多邊形，也可以將

38、它們組合使用。(2)boundary模式：在各個(gè)邊界段上給出邊界條件。(3)pde模式：確定方程的類型、c,a,f和d.也能夠在不同的子區(qū)域上設(shè)置不同的系數(shù)(反映材料的性質(zhì)).用gui設(shè)置定解問題包括三個(gè)模式(mode):(1)mesh模式：生成網(wǎng)格、自動(dòng)控制參數(shù)。 (2)solve模式：對(duì)于橢圓型方程還能求非線性和自適應(yīng)解。對(duì)于拋物型和雙曲型方程，設(shè)置初始邊值條件后能求出給定t時(shí)刻的解。求解后可以加密網(wǎng)格再求解。為此，選擇pde菜單中pde mode命令，進(jìn)入pde模式(如圖12)，此時(shí)需要確定方程的類型、c,a,f和d.也能夠在不同的子區(qū)域上設(shè)置不同的系數(shù)(反映材料的性質(zhì))：區(qū)域4鍵入mu

39、=1,j=1;區(qū)域1、2均鍵入這里系數(shù)含有非線性；區(qū)域3鍵入mu=1,j=0.圖12 (5)如何初始化網(wǎng)格： 選擇mesh菜單中initialize mesh命令，即可對(duì)平面圖形進(jìn)行網(wǎng)格剖分. (6)如何進(jìn)行網(wǎng)格加密： 選擇mesh菜單中refine mesh命令，即可實(shí)現(xiàn)對(duì)網(wǎng)格進(jìn)行加密. (7)怎樣進(jìn)行求解，如何顯示圖形： 選擇solve菜單中solve pde命令，解偏微分方程并顯示圖形解(如圖13所示),圖13(8)顯示三維圖形解： 單擊plot菜單中parameters.選項(xiàng)，打開plot selection對(duì)話框，選中color,contour ,height(3-d plot)和a

40、rrows四項(xiàng)，設(shè)置contour plot levels:10,然后單擊plot按鈕，顯示三維圖形解.如圖14所示圖14 解的3維圖形(2)拋物型方程：考慮拋物線偏微分方程：它描述平面的熱傳輸現(xiàn)象，以及由軸對(duì)稱三維問題經(jīng)過降維后的熱傳輸問題，其中t為溫度，其他參數(shù)為： 在pde toolbox中我們可以做溫度、溫度梯度、熱流的可視化圖形。 考慮一個(gè)方形區(qū)域，導(dǎo)熱系數(shù)為10，密度為2。在方形區(qū)域內(nèi)有一個(gè)圓形區(qū)域的熱源為4，導(dǎo)熱系數(shù)為2，密度為1。兩個(gè)區(qū)域的比熱都為0.1。以下我們?cè)敿?xì)說明利用gui求解的過程，依次回答以下問題：(1)如何啟動(dòng)pde toolbox工具箱：?jiǎn)?dòng)matlab7.0，

41、在其命令窗口(command window)中鍵入命令pdetool，這時(shí)系統(tǒng)將立即產(chǎn)生偏微分方程工具箱(pde toolbox)的圖形用戶界面(graphical user interface,簡(jiǎn)記為gui)，也就實(shí)現(xiàn)了pde toolbox工具箱的啟動(dòng)。此外，由于此問題所涉及的是平面熱傳輸?shù)膯栴}，因此，我們需要將系統(tǒng)默認(rèn)的應(yīng)用模型模式選定為heat transfer. (2)怎樣分步完成平面幾何造型： 為了便于平面圖形的觀察，為此，我們對(duì)坐標(biāo)的范圍顯示進(jìn)行設(shè)置，即單擊options菜單中的grid spacing對(duì)話框中選auto;打開axes limits對(duì)話框，設(shè)置x軸和y軸的范圍為0,3,單擊apply即可完成設(shè)置。 此外,為了使兩坐標(biāo)軸的尺寸比例一致,我們重新調(diào)整坐標(biāo)軸的范圍，即單擊options菜單中axes equal選項(xiàng)。 在pde toolbox面板中依次拖出圓方形分別雙擊所繪制的圖形，對(duì)于圓:設(shè)置x-center為1.5，y-center為1.5，radius為1;對(duì)于方形;設(shè)置left為0，bottom分別為0,width分別為3，height分別為3. 這樣便可以得到圖15所示的平面圖.圖15 因?yàn)樗型?/p>