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

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

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

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

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

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

6、件廣義諾依曼（Neumann）邊界條件齊次邊界條件： g=0，r =0第一類邊界條件： Diriclet邊界條件第三類邊界條件： Neumann邊界條件第二類邊界條件： q=0偏微分方程工具箱求解定解問題邊界條件：Page 8狄利克里（偏微分方程工具箱求解定解問題第四步：設置方程類型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第八步：繪制等值線圖和矢量場圖 Plot Parameter，選擇Contour和ArrowsPlotPage 9偏微分方程工具箱求解定解問題第四步：設置方程類型Page 9雙曲型：波動方程的求解與可視化雙曲型：波動方程的求解與靜態(tài)（或動態(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)格坐標描述矩陣決定方程的類型時間矩陣邊界條件雙曲型：波動方程的求解與可

8、視化雙曲型：波動方程的求解與靜態(tài)（雙曲型：波動方程的求解與可視化網(wǎng)格初始化命令：(1) p,e,t=initmesh(g)將求解區(qū)域進行三角形網(wǎng)格化，輸出的p、e、t是網(wǎng)格數(shù)據(jù)p-描述網(wǎng)格中點的x、y坐標e-邊緣矩陣，t-三角矩陣，描述區(qū)域的頂點g-描述求解區(qū)域幾何形狀(2) p,e,t=refinemesh(g,p,e,t) 迭代過程，得到更細小的網(wǎng)格，使結(jié)果更精確Page 11雙曲型：波動方程的求解與可視化網(wǎng)格初始化命令：Page 11雙曲型：波動方程的求解與可視化2.動畫圖形顯示 將所得的解形象地表示出來，為了加速繪圖，首先把三角形網(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) 用此命令之前，應先用一個tri2grid命令得出矩陣tn、a2、a3用此方法可以加快速度Page 12三角形網(wǎng)格的矩陣矩形網(wǎng)格的坐標點各時刻三角形網(wǎng)格中的解插值法求得矩形網(wǎng)格點上的u值內(nèi)插法的系數(shù)格點的指針矩陣雙曲型：波動方程的求解與可視化2.動畫圖形顯示 Page 雙曲型：波動方程的求解與可視化主要的繪圖（包括動畫）命令函數(shù)有：moviein、movie、pedplot、pdesurf等Page 13雙曲型：波動方程的

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

11、reg); （3）初始條件 x=p(1,:); % 注意坐標向量都是列向量 y=p(2,:); u0=atan(sin(pi/2*x); ut0=2*cos(pi*x).*exp(cos(pi/2*y);Page 15雙曲型：波動方程的求解與可視化【解】采用步驟如下Page 1雙曲型：波動方程的求解與可視化（4）在時間段05內(nèi)的31個點上求解 n=31; tlist=linspace(0,5,n); % 在05之間產(chǎn)生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雙曲型：波動方程的求解與可視化（4）在時間段05內(nèi)的31個雙曲型：波動方程的求解與可視化 解u1的動態(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時間內(nèi)動畫顯示 newplot; %建立新的坐標系 M=moviein(n); umax=max(max(u1); umin=min(min(u1);Page 17雙曲型：波動方程的求解與可視化 解u1的動態(tài)圖形顯示雙曲型：波動方程的求解與可視化for i=1: nif rem(i,10) = 0 %當n是10的整數(shù)倍時， fprintf(%d,i); %在命令窗口打印出相應的數(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雙曲型：波動方程的求解與可視化for i=1: nPage雙曲型：波動方程的求解與可視化%運行結(jié)果是：10 20 30 done動態(tài)解圖可以直接通過MATLAB仿真程序執(zhí)行看出，圖1是動態(tài)圖的某一瞬間的解的分布。P

15、age 19雙曲型：波動方程的求解與可視化%運行結(jié)果是：Page 19雙曲型：波動方程的求解與可視化Page 20圖1 某一瞬時的波動方程的解圖雙曲型：波動方程的求解與可視化Page 20圖1 某一瞬時的雙曲型：波動方程的求解與可視化Page 21例2 討論弦的一端x=0,固定，x=L一端受迫作諧振動2sint，弦的初始位移和初始速度為零，給出弦振動的解圖?！窘狻?根據(jù)題意得定解問題為： 解析解為： 雙曲型：波動方程的求解與可視化Page 21例2 討論弦的一雙曲型：波動方程的求解與可視化計算機仿真程序如下：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雙曲型：波動方程的求解與可視化計算機仿真程序如下：Page 雙曲型：波動方程的求解與可視化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雙曲型：波動方程的求解與可視化figure(1)Page 2雙曲型：波動方程的求解與可視化Page 24圖2 波動方程解析解的分布雙曲型：波動方程的求解與可視化Page 24圖2 波動方程解拋物型：熱傳導方程的求解與可視化熱傳導方程屬于拋物線方程，在MATLAB中是指如下形式：求解拋物線方

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

20、 er; fprintf(Error: %e. Number of nodes: %dn,er,size(p,2); end % 運行結(jié)果是 Error: 1.292265e-002. Number of nodes: 25 Page 37橢圓型：穩(wěn)定場方程的求解與可視化(3) 迭代直至得到誤差允許橢圓型：穩(wěn)定場方程的求解與可視化(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)定場方程的求解與可視化(4) 把結(jié)果用圖形表示：P橢圓型：穩(wěn)定場方程的求解與可視化Page 39圖5 精確解與仿真解的誤差解圖(a) 泊松方程的解析解圖(b)精確解與仿真解的誤差解橢圓型：穩(wěn)定場方程的求解與可視化Page 39圖5 精確解與橢圓型：穩(wěn)定場方程的求解與可視化 例6 在矩形區(qū)域0 xa，-b/2yb/2上，對滿足泊松方程 ，且邊界上的值為零的定解問題的解，給出計算機仿真圖形。 【解】所討論的定解問題即為：Page 40橢圓型：穩(wěn)定場方程的求解與可視化 例6 在矩形區(qū)域0 xa橢圓型

22、：穩(wěn)定場方程的求解與可視化 解析解的表達式： 計算機仿真程序：（程序中取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)定場方程的求解與可視化 解析解的表達式：Pa橢圓型：穩(wěn)定場方程的求解與可視化Z=Z1+X.*Y.*(a3-X.3)/12; colorma

23、p(hot); mesh(X,Y,Z) view(119,7);Page 42圖6 矩形域的泊松方程解的分布橢圓型：穩(wěn)定場方程的求解與可視化Z=Z1+X.*Y.*(aPage 43本征值問題的求解與可視化二維本征值問題矩形區(qū)域的本征模與本征振動邊長為b和c的的四周固定的矩形膜振動的本征值問題為采用分離變量法可以得到本征模和本征值為Page 43本征值問題的求解與可視化二維本征值問題邊長為bPage 44本征值問題的求解與可視化繪制前4個本征函數(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個本征函數(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個本征函數(shù)的圖本征值問題的求解與可視化Page 45圖7 前4個本征函數(shù)的圖形本征值問題的求解與可視化Page 45圖7 前4個本征函數(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