數(shù)學(xué)物理方程的計(jì)算機(jī)求解和可視化課件

1、數(shù)學(xué)物理建模與計(jì)算機(jī)輔助設(shè)計(jì)專題4：數(shù)學(xué)物理方程的計(jì)算機(jī)求解和可視化數(shù)學(xué)物理建模與計(jì)算機(jī)輔助設(shè)計(jì)專題4：數(shù)學(xué)物理方程的計(jì)算機(jī)求解Page 2本專題主要內(nèi)容與參考資料 主要內(nèi)容偏微分方程的計(jì)算機(jī)仿真求解方法雙曲型（Hyperbolic）：波動(dòng)方程的求解與可視化拋物型（Parabolic）：熱傳導(dǎo)方程的求解與可視化橢圓型（Elliptic）：穩(wěn)定場(chǎng)方程的求解與可視化特征值問題的求解與可視化利用Matlab的PDE工具箱求解并進(jìn)行可視化 參考資料楊華軍, 數(shù)學(xué)物理方法, 電子工業(yè)出版社彭芳麟, 數(shù)學(xué)物理方程的Matlab解法與可視化, 清華大學(xué)出版社Page 2本專題主要內(nèi)容與參考資料 主要內(nèi)容仿真

2、求解偏微分方程的類型通用的線性二階偏微分方程：偏微分方程類型分為：雙曲型方程： 拋物型方程： 橢圓型方程： 特征值問題： 特征值偏微分方程中不含參數(shù)f .Page 3仿真求解偏微分方程的類型通用的線性二階偏微分方程：Page 偏微分方程的仿真求解方法偏微分方程的計(jì)算機(jī)仿真求解方法：（1）MATLAB的偏微分方程工具箱（PDE Toolbox） 有限元法（2）MATLAB仿真，M文件編程 典型偏微分的解的靜態(tài)（或動(dòng)態(tài)）三維可視化。Page 4偏微分方程的仿真求解方法偏微分方程的計(jì)算機(jī)仿真求解方法：Pa有限元法定義將連續(xù)的求解域離散成一組有限個(gè)、按照一定方式相互連結(jié)在一起的單元的組合體。將PDE轉(zhuǎn)

3、換成離散的線性代數(shù)方程系統(tǒng)進(jìn)行求解。特點(diǎn)各種復(fù)雜單元可以用來(lái)對(duì)復(fù)雜的幾何形狀的求解域進(jìn)行模型化處理。各節(jié)點(diǎn)上的解的近似函數(shù)可以用來(lái)求解整個(gè)求解域上任意點(diǎn)的結(jié)果。Page 5有限元法定義Page 5MATLAB的偏微分方程工具箱（PDE Toolbox）用圖形用戶界面（Graphical User Interface，簡(jiǎn)記作GUI）求解PDE問題主要采用以下三個(gè)步驟：設(shè)置PDE的定解問題 即設(shè)置二維定解區(qū)域、邊界條件以及方程的形式和系數(shù)；(2) 用有限元法（FEM）求解PDE即網(wǎng)格的生成、方程的離散以及求出數(shù)值解； Mesh：生成網(wǎng)格，自動(dòng)控制網(wǎng)格參數(shù) Solve：求解設(shè)置初始邊值條件后，能給出

4、t時(shí)刻的解；可以求出區(qū)間內(nèi)的特征值求解后可以加密網(wǎng)格再求解(3) 解的可視化從GUI使用Plot方法實(shí)現(xiàn)可視化用Color、Height、Vector等作圖對(duì)于含時(shí)方程，還可以生成解的動(dòng)畫用PDE Toolbox可以求解的基本方程有：橢圓方程、拋物方程、雙曲方程、特征值方程、橢圓方程組以及非線性橢圓方程。Page 6MATLAB的偏微分方程工具箱（PDE Toolbox）用圖偏微分方程工具箱求解定解問題例1 解熱傳導(dǎo)方程 邊界條件是齊次類型（ u=0），定解區(qū)域自定?！窘狻康谝徊剑?jiǎn)?dòng)MATLAB，鍵入命令pdetool并回車，就進(jìn)入GUI Options Grid，打開柵格第二步：選定定解區(qū)

5、域 繪制橢圓E1、圓E2、矩形R1、圓E3; 在Set formula欄中鍵入E1-E2+R1-E3.第三步：選取邊界Boundary Boundary Mode，進(jìn)入邊界模式；Boundary Remove All Subdomain Borders，去掉子域邊界;Boundary Specify Boundary Conditions Boundary Conditions 齊次Diriclet條件.Page 7偏微分方程工具箱求解定解問題例1 解熱傳導(dǎo)方程 Pag偏微分方程工具箱求解定解問題邊界條件：解方程所需要的邊界條件可以是以下兩種形式：Page 8狄利克里（Diriclet）邊界條

6、件廣義諾依曼（Neumann）邊界條件齊次邊界條件： g=0，r =0第一類邊界條件： Diriclet邊界條件第三類邊界條件： Neumann邊界條件第二類邊界條件： q=0偏微分方程工具箱求解定解問題邊界條件：Page 8狄利克里（偏微分方程工具箱求解定解問題第四步：設(shè)置方程類型PDE ModePDE Secification，選擇方程類型:拋物型第五步：劃分網(wǎng)格 Mesh Initialize Mesh，網(wǎng)格剖分； Mesh Refine Mesh，網(wǎng)格密集化.第六步： 解偏微分方程并顯示圖形解 Solve Solve PDE，解偏微分方程并顯示圖形解。第七步：三維可視化 Plot Pa

7、rameter，選擇Color, Height(3-D plot)，Show mesh第八步：繪制等值線圖和矢量場(chǎng)圖 Plot Parameter，選擇Contour和ArrowsPlotPage 9偏微分方程工具箱求解定解問題第四步：設(shè)置方程類型Page 9雙曲型：波動(dòng)方程的求解與可視化雙曲型：波動(dòng)方程的求解與靜態(tài)（或動(dòng)態(tài)）三維可視化 1. 求解雙曲型方程求解雙曲型方程調(diào)用形式如下：u1=hyperbolic（u0,ut0,tlist,b,p,e,t,c,a,f,d）Page 10（a、c、d、f是參數(shù)）初始條件ut即是網(wǎng)格坐標(biāo)描述矩陣決定方程的類型時(shí)間矩陣邊界條件雙曲型：波動(dòng)方程的求解與可

8、視化雙曲型：波動(dòng)方程的求解與靜態(tài)（雙曲型：波動(dòng)方程的求解與可視化網(wǎng)格初始化命令：(1) p,e,t=initmesh(g)將求解區(qū)域進(jìn)行三角形網(wǎng)格化，輸出的p、e、t是網(wǎng)格數(shù)據(jù)p-描述網(wǎng)格中點(diǎn)的x、y坐標(biāo)e-邊緣矩陣，t-三角矩陣，描述區(qū)域的頂點(diǎn)g-描述求解區(qū)域幾何形狀(2) p,e,t=refinemesh(g,p,e,t) 迭代過(guò)程，得到更細(xì)小的網(wǎng)格，使結(jié)果更精確Page 11雙曲型：波動(dòng)方程的求解與可視化網(wǎng)格初始化命令：Page 11雙曲型：波動(dòng)方程的求解與可視化2.動(dòng)畫圖形顯示 將所得的解形象地表示出來(lái)，為了加速繪圖，首先把三角形網(wǎng)格轉(zhuǎn)化成矩形網(wǎng)格調(diào)用形式如下： (1) uxy=tri

9、2grid（p，t，u1，x，y） (2) uxy,tn,a2,a3=tri2grid(p,t,u,x,y) (3) uxy=tri2grid(p,t,u,tn,a2,a3) 用此命令之前，應(yīng)先用一個(gè)tri2grid命令得出矩陣tn、a2、a3用此方法可以加快速度Page 12三角形網(wǎng)格的矩陣矩形網(wǎng)格的坐標(biāo)點(diǎn)各時(shí)刻三角形網(wǎng)格中的解插值法求得矩形網(wǎng)格點(diǎn)上的u值內(nèi)插法的系數(shù)格點(diǎn)的指針矩陣雙曲型：波動(dòng)方程的求解與可視化2.動(dòng)畫圖形顯示 Page 雙曲型：波動(dòng)方程的求解與可視化主要的繪圖（包括動(dòng)畫）命令函數(shù)有：moviein、movie、pedplot、pdesurf等Page 13雙曲型：波動(dòng)方程的

10、求解與可視化主要的繪圖（包括動(dòng)畫）命令函數(shù)雙曲型：波動(dòng)方程的求解與可視化例1 用MATLAB求解下面波動(dòng)方程定解問題并動(dòng)態(tài)顯示解的分布 方法1：其解可以用偏微分方程工具箱求得，繪制出其圖形。 方法2：求出解析解，再利用MATLAB編程繪制出其解析的圖形分布。Page 14雙曲型：波動(dòng)方程的求解與可視化例1 用MATLAB求解下面波雙曲型：波動(dòng)方程的求解與可視化【解】采用步驟如下（1）題目定義 g=squareg; % 定義單位方形區(qū)域 b=squareb3; % 左右零邊界條件，頂?shù)琢銓?dǎo)數(shù)邊界條件 c=1;a=0;f=0;d=1; （2）初始的粗糙網(wǎng)格化 p,e,t=initmesh(squa

11、reg); （3）初始條件 x=p(1,:); % 注意坐標(biāo)向量都是列向量 y=p(2,:); u0=atan(sin(pi/2*x); ut0=2*cos(pi*x).*exp(cos(pi/2*y);Page 15雙曲型：波動(dòng)方程的求解與可視化【解】采用步驟如下Page 1雙曲型：波動(dòng)方程的求解與可視化（4）在時(shí)間段05內(nèi)的31個(gè)點(diǎn)上求解 n=31; tlist=linspace(0,5,n); % 在05之間產(chǎn)生n個(gè)均勻的時(shí)間點(diǎn)（5）求解此雙曲問題u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d); 得到如下結(jié)果：Page 16%428 success

12、ful steps%62 failed attempts%982 function evaluations%1 partial derivatives%142 LU decompositions%981 solutions of linear systems%Time: 0.166667%Time: 0.333333%Time: 0.5%Time: 4.5%Time: 4.66667%Time: 4.83333%Time: 5雙曲型：波動(dòng)方程的求解與可視化（4）在時(shí)間段05內(nèi)的31個(gè)雙曲型：波動(dòng)方程的求解與可視化 解u1的動(dòng)態(tài)圖形顯示：（6）矩形網(wǎng)格插值 delta=-1:0.1:1; uxy

13、, tn, a2, a3=tri2grid(p, t, u1(:, 1), delta, delta); gp=tn; a2; a3;（7）在05時(shí)間內(nèi)動(dòng)畫顯示 newplot; %建立新的坐標(biāo)系 M=moviein(n); umax=max(max(u1); umin=min(min(u1);Page 17雙曲型：波動(dòng)方程的求解與可視化 解u1的動(dòng)態(tài)圖形顯示雙曲型：波動(dòng)方程的求解與可視化for i=1: nif rem(i,10) = 0 %當(dāng)n是10的整數(shù)倍時(shí)， fprintf(%d,i); %在命令窗口打印出相應(yīng)的數(shù)字endpdeplot(p, e, t, xydata, u1(:, i

14、),zdata, u1(:,i), zstyle, continuous, mesh, on,xygrid,on, gridparam, gp, colorbar, off); axis(-1, 1, -1, 1 umin umax); caxis(umin umax); M(:, i)=getframe; if i =n fprintf(donen ); endendPage 18雙曲型：波動(dòng)方程的求解與可視化for i=1: nPage雙曲型：波動(dòng)方程的求解與可視化%運(yùn)行結(jié)果是：10 20 30 done動(dòng)態(tài)解圖可以直接通過(guò)MATLAB仿真程序執(zhí)行看出，圖1是動(dòng)態(tài)圖的某一瞬間的解的分布。P

15、age 19雙曲型：波動(dòng)方程的求解與可視化%運(yùn)行結(jié)果是：Page 19雙曲型：波動(dòng)方程的求解與可視化Page 20圖1 某一瞬時(shí)的波動(dòng)方程的解圖雙曲型：波動(dòng)方程的求解與可視化Page 20圖1 某一瞬時(shí)的雙曲型：波動(dòng)方程的求解與可視化Page 21例2 討論弦的一端x=0,固定，x=L一端受迫作諧振動(dòng)2sint，弦的初始位移和初始速度為零，給出弦振動(dòng)的解圖?！窘狻?根據(jù)題意得定解問題為： 解析解為： 雙曲型：波動(dòng)方程的求解與可視化Page 21例2 討論弦的一雙曲型：波動(dòng)方程的求解與可視化計(jì)算機(jī)仿真程序如下：clear a=1;l=1;A=2.0;w=6;x=0:0.05:1;t=0:0.00

16、1:4.3;X,T=meshgrid(x,t);u0=A*sin(w*X./a).*sin(w.*T)/sin(w*l/a);u=0;for n=1:100; uu=(-1)(n+1)*sin(n*pi*X/l).*sin(n*pi*a*T/l)/(w*w/a/a- n*n*pi*pi/l/l); u=u+uu;endu=u0+2*A*w/a/l.*u;Page 22雙曲型：波動(dòng)方程的求解與可視化計(jì)算機(jī)仿真程序如下：Page 雙曲型：波動(dòng)方程的求解與可視化figure(1)axis(0 1 -5 5)h=plot(x,u(1,:),linewidth,3);set(h,erasemode,xo

17、r);for j=2:length(t); set(h,ydata,u(j,:); axis(0 1 -5 5) drawnowendfigure(2)waterfall(X(1:50:3000,:),T(1:50:3000,:),u(1:50:3000,:)Xlabel(x)Ylabel(t)Page 23雙曲型：波動(dòng)方程的求解與可視化figure(1)Page 2雙曲型：波動(dòng)方程的求解與可視化Page 24圖2 波動(dòng)方程解析解的分布雙曲型：波動(dòng)方程的求解與可視化Page 24圖2 波動(dòng)方程解拋物型：熱傳導(dǎo)方程的求解與可視化熱傳導(dǎo)方程屬于拋物線方程，在MATLAB中是指如下形式：求解拋物線方

18、程使用如下命令：u=parabolic(u0,ut0, tlist, b, p, e, t, c, a, f, d) parabolic函數(shù)性質(zhì)與hyperbolic大致相同Page 25初始條件ut即是網(wǎng)格坐標(biāo)描述矩陣決定方程的類型時(shí)間矩陣邊界條件拋物型：熱傳導(dǎo)方程的求解與可視化熱傳導(dǎo)方程屬于拋物線方程，在拋物型：熱傳導(dǎo)方程的求解與可視化例3 求解下列熱傳導(dǎo)定解問題求解域是方形區(qū)域，其中空間坐標(biāo)的個(gè)數(shù)由具體問題確定Page 26拋物型：熱傳導(dǎo)方程的求解與可視化例3 求解下列熱傳導(dǎo)定解問題拋物型：熱傳導(dǎo)方程的求解與可視化【解】步驟如下：(1) 題目定義 g=squareg; % 定義單位方形區(qū)

19、域 b=squareb1; % 定義零邊界條件 c=1;a=0;f=1;d=1;(2) 網(wǎng)格化 p,e,t=initmesh(g);(3) 定義初始條件 u0=zeros(size(p,2),1); %產(chǎn)生零矩陣u0，size(p,2)返回p的列數(shù) ix=find(sqrt(p(1,:).2+p(2,:).2)0.001, p,e,t=refinemesh(g,p,e,t); u=assempde(b,p,e,t,c,a,f); exact=(1-p(1,:).2-p(2,:).2)/4; er=norm(u-exact,inf); % er是u-exact的無(wú)窮大的模 error=error

20、 er; fprintf(Error: %e. Number of nodes: %dn,er,size(p,2); end % 運(yùn)行結(jié)果是 Error: 1.292265e-002. Number of nodes: 25 Page 37橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化(3) 迭代直至得到誤差允許橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化(4) 把結(jié)果用圖形表示：%pdesurf(p,t,u); pdeplot (p,e,t,xydata, u, zdata, u, mesh, on); figure; %pdesurf(p,t,u-exact);pdeplot (p,e,t,xydata,u-exa

21、ct, zdata,u-exact, mesh, on); % 誤差解圖顯示 Page 38橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化(4) 把結(jié)果用圖形表示：P橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化Page 39圖5 精確解與仿真解的誤差解圖(a) 泊松方程的解析解圖(b)精確解與仿真解的誤差解橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化Page 39圖5 精確解與橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化 例6 在矩形區(qū)域0 xa，-b/2yb/2上，對(duì)滿足泊松方程 ，且邊界上的值為零的定解問題的解，給出計(jì)算機(jī)仿真圖形。 【解】所討論的定解問題即為：Page 40橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化 例6 在矩形區(qū)域0 xa橢圓型

22、：穩(wěn)定場(chǎng)方程的求解與可視化 解析解的表達(dá)式： 計(jì)算機(jī)仿真程序：（程序中取a=8，b=8 ）syms a ba=8; b=8;X,Y=meshgrid(0:0.2:a,-b/2:0.2:b/2)Z1=0;for n=1:1:10 Z2=a4*b*(-1)n*n2*pi2+2-2*(-1)n)*sinh(n*pi.*Y/a).* sin(n*pi.*X/a)/(n5*pi5*sinh(n*pi*b/(2*a); Z1=Z1+Z2;end Page 41橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化 解析解的表達(dá)式：Pa橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化Z=Z1+X.*Y.*(a3-X.3)/12; colorma

23、p(hot); mesh(X,Y,Z) view(119,7);Page 42圖6 矩形域的泊松方程解的分布橢圓型：穩(wěn)定場(chǎng)方程的求解與可視化Z=Z1+X.*Y.*(aPage 43本征值問題的求解與可視化二維本征值問題矩形區(qū)域的本征模與本征振動(dòng)邊長(zhǎng)為b和c的的四周固定的矩形膜振動(dòng)的本征值問題為采用分離變量法可以得到本征模和本征值為Page 43本征值問題的求解與可視化二維本征值問題邊長(zhǎng)為bPage 44本征值問題的求解與可視化繪制前4個(gè)本征函數(shù)的圖形%P70_1.ma=2; b=1;m,n=meshgrid(1:3);L=(m*pi./b).2+(n*pi./b).2) %求本征值x=0:0.

24、01:a; y=0:0.01:b;X,Y=meshgrid(x,y);w11=sin(pi*Y./b).*sin(pi*X./a); w12=sin(2*pi*Y./b).*sin(pi*X./a); w21=sin(pi*Y./b).*sin(2*pi*X./a);w22=sin(pi*Y./b).*sin(3*pi*X./a);%求前4個(gè)本征函數(shù)figuresubplot(2,2,1); mesh(X,Y,w11); subplot(2,2,2); mesh(X,Y,w12);subplot(2,2,3); mesh(X,Y,w21); subplot(2,2,4); mesh(X,Y,w

25、22);L = 19.7392 49.3480 98.6960 49.3480 78.9568 128.3049 98.6960 128.3049 177.6529Page 44本征值問題的求解與可視化繪制前4個(gè)本征函數(shù)的圖本征值問題的求解與可視化Page 45圖7 前4個(gè)本征函數(shù)的圖形本征值問題的求解與可視化Page 45圖7 前4個(gè)本征函數(shù)的本征值問題的求解與可視化figure(2)subplot(2,2,1); mesh(X,Y,w11); view(0,90)subplot(2,2,2); mesh(X,Y,w12); view(0,90)subplot(2,2,3); mesh(X,

26、Y,w21); view(0,90)subplot(2,2,4); mesh(X,Y,w22); view(0,90)Page 46本征值問題的求解與可視化figure(2)Page 46本征值問題的求解與可視化figure(3)subplot(2,2,1); contour(X,Y,w11); subplot(2,2,2); contour(X,Y,w12); subplot(2,2,3); contour(X,Y,w21); subplot(2,2,4); contour(X,Y,w22); Page 47本征值問題的求解與可視化figure(3)Page 47本征值問題的求解與可視化比較

27、下面兩副圖的區(qū)別：figure(4)subplot(2,2,1); C,h=contour(X,Y,w11,20); clabel(C,h,manual);subplot(2,2,2); C,h=contour(X,Y,w12,20); clabel(C,h,manual); subplot(2,2,3); C,h=contour(X,Y,w21,20); clabel(C,h,manual); subplot(2,2,4); C,h=contour(X,Y,w22,20); clabel(C,h,manual); figure(5)subplot(2,2,1); C,h=contour(X,Y