課堂授課專題4：數(shù)學(xué)物理方程的計(jì)算機(jī)求解和可視化

數(shù)學(xué)物理建模與計(jì)算機(jī)輔助設(shè)計(jì)專題4：數(shù)學(xué)物理方程的計(jì)算機(jī)求解和可視化Page2本專題主要內(nèi)容與參考資料

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

拋物型方程：

橢圓型方程：

特征值問題：

特征值偏微分方程中不含參數(shù)f.Page3偏微分方程的仿真求解方法偏微分方程的計(jì)算機(jī)仿真求解方法：（1）MATLAB的偏微分方程工具箱（PDEToolbox）

有限元法（2）MATLAB仿真，M文件編程

典型偏微分的解的靜態(tài)（或動(dòng)態(tài)）三維可視化。Page4有限元法定義將連續(xù)的求解域離散成一組有限個(gè)、按照一定方式相互連結(jié)在一起的單元的組合體。將PDE轉(zhuǎn)換成離散的線性代數(shù)方程系統(tǒng)進(jìn)行求解。特點(diǎn)各種復(fù)雜單元可以用來對復(fù)雜的幾何形狀的求解域進(jìn)行模型化處理。各節(jié)點(diǎn)上的解的近似函數(shù)可以用來求解整個(gè)求解域上任意點(diǎn)的結(jié)果。Page5MATLAB的偏微分方程工具箱（PDEToolbox）用圖形用戶界面（GraphicalUserInterface，簡記作GUI）求解PDE問題主要采用以下三個(gè)步驟：設(shè)置PDE的定解問題．

即設(shè)置二維定解區(qū)域、邊界條件以及方程的形式和系數(shù)；(2)用有限元法（FEM）求解PDE．即網(wǎng)格的生成、方程的離散以及求出數(shù)值解；

Mesh：生成網(wǎng)格，自動(dòng)控制網(wǎng)格參數(shù)．Solve：求解．設(shè)置初始邊值條件后，能給出t時(shí)刻的解；可以求出區(qū)間內(nèi)的特征值．求解后可以加密網(wǎng)格再求解．(3)解的可視化．從GUI使用Plot方法實(shí)現(xiàn)可視化．用Color、Height、Vector等作圖．對于含時(shí)方程，還可以生成解的動(dòng)畫．用PDEToolbox可以求解的基本方程有：橢圓方程、拋物方程、雙曲方程、特征值方程、橢圓方程組以及非線性橢圓方程。．Page6偏微分方程工具箱求解定解問題例1解熱傳導(dǎo)方程

邊界條件是齊次類型（

u=0），定解區(qū)域自定。【解】第一步：啟動(dòng)MATLAB，鍵入命令pdetool并回車，就進(jìn)入GUI．

Options→Grid，打開柵格．第二步：選定定解區(qū)域繪制橢圓E1、圓E2、矩形R1、圓E3;在Setformula欄中鍵入E1-E2+R1-E3.第三步：選取邊界Boundary→BoundaryMode，進(jìn)入邊界模式；Boundary→RemoveAllSubdomainBorders，去掉子域邊界;Boundary→SpecifyBoundaryConditions

→BoundaryConditions→齊次Diriclet條件.Page7偏微分方程工具箱求解定解問題邊界條件：解方程所需要的邊界條件可以是以下兩種形式：Page8狄利克里（Diriclet）邊界條件廣義諾依曼（Neumann）邊界條件齊次邊界條件：g=0，r=0第一類邊界條件：Diriclet邊界條件第三類邊界條件：Neumann邊界條件第二類邊界條件：q=0偏微分方程工具箱求解定解問題第四步：設(shè)置方程類型PDEMode→PDESecification，選擇方程類型:拋物型第五步：劃分網(wǎng)格

Mesh→InitializeMesh，網(wǎng)格剖分；

Mesh→RefineMesh，網(wǎng)格密集化.第六步：解偏微分方程并顯示圖形解

Solve→SolvePDE，解偏微分方程并顯示圖形解。第七步：三維可視化

Plot

→Parameter，選擇Color,Height(3-Dplot)，Showmesh第八步：繪制等值線圖和矢量場圖

Plot→Parameter，選擇Contour和Arrows．PlotPage9雙曲型：波動(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）Page10（a、c、d、f是參數(shù)）初始條件ut即是網(wǎng)格坐標(biāo)描述矩陣決定方程的類型時(shí)間矩陣邊界條件雙曲型：波動(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)

迭代過程，得到更細(xì)小的網(wǎng)格，使結(jié)果更精確．Page11雙曲型：波動(dòng)方程的求解與可視化2.動(dòng)畫圖形顯示

將所得的解形象地表示出來，為了加速繪圖，首先把三角形網(wǎng)格轉(zhuǎn)化成矩形網(wǎng)格．調(diào)用形式如下：

(1)uxy=tri2grid（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．用此方法可以加快速度．Page12三角形網(wǎng)格的矩陣矩形網(wǎng)格的坐標(biāo)點(diǎn)各時(shí)刻三角形網(wǎng)格中的解插值法求得矩形網(wǎng)格點(diǎn)上的u值內(nèi)插法的系數(shù)格點(diǎn)的指針矩陣雙曲型：波動(dòng)方程的求解與可視化主要的繪圖（包括動(dòng)畫）命令函數(shù)有：moviein、movie、pedplot、pdesurf等Page13雙曲型：波動(dòng)方程的求解與可視化例1用MATLAB求解下面波動(dòng)方程定解問題并動(dòng)態(tài)顯示解的分布

方法1：其解可以用偏微分方程工具箱求得，繪制出其圖形。

方法2：求出解析解，再利用MATLAB編程繪制出其解析的圖形分布。Page14雙曲型：波動(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('squareg');（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));Page15雙曲型：波動(dòng)方程的求解與可視化（4）在時(shí)間段0～5內(nèi)的31個(gè)點(diǎn)上求解n=31;

tlist=linspace(0,5,n);%在0～5之間產(chǎn)生n個(gè)均勻的時(shí)間點(diǎn)（5）求解此雙曲問題u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);得到如下結(jié)果：Page16%428successfulsteps%62failedattempts%982functionevaluations%1partialderivatives%142LUdecompositions%981solutionsoflinearsystems%Time:0.166667%Time:0.333333%Time:0.5%…%Time:4.5%Time:4.66667%Time:4.83333%Time:5雙曲型：波動(dòng)方程的求解與可視化解u1的動(dòng)態(tài)圖形顯示：（6）矩形網(wǎng)格插值

delta=-1:0.1:1;[uxy,tn,a2,a3]=tri2grid(p,t,u1(:,1),delta,delta);gp=[tn;a2;a3];（7）在0～5時(shí)間內(nèi)動(dòng)畫顯示newplot;%建立新的坐標(biāo)系M=moviein(n);umax=max(max(u1));umin=min(min(u1));Page17雙曲型：波動(dòng)方程的求解與可視化fori=1:nifrem(i,10)==0%當(dāng)n是10的整數(shù)倍時(shí)，fprintf('%d’,i);%在命令窗口打印出相應(yīng)的數(shù)字endpdeplot(p,e,t,'xydata',u1(:,i),'zdata',u1(:,i),'zstyle','continuous','mesh','on','xygrid','on','gridparam',gp,'colorbar','off');axis([-1,1,-1,1uminumax]);caxis([uminumax]);M(:,i)=getframe;ifi==nfprintf('done\n');endendPage18雙曲型：波動(dòng)方程的求解與可視化%運(yùn)行結(jié)果是：102030done動(dòng)態(tài)解圖可以直接通過MATLAB仿真程序執(zhí)行看出，圖1是動(dòng)態(tài)圖的某一瞬間的解的分布。Page19雙曲型：波動(dòng)方程的求解與可視化Page20圖1某一瞬時(shí)的波動(dòng)方程的解圖雙曲型：波動(dòng)方程的求解與可視化Page21例2討論弦的一端x=0,固定，x=L一端受迫作諧振動(dòng)2sinωt，弦的初始位移和初始速度為零，給出弦振動(dòng)的解圖。【解】

根據(jù)題意得定解問題為：

解析解為：

雙曲型：波動(dòng)方程的求解與可視化計(jì)算機(jī)仿真程序如下：cleara=1;l=1;A=2.0;w=6;x=0:0.05:1;t=0:0.001:4.3;[X,T]=meshgrid(x,t);u0=A*sin(w*X./a).*sin(w.*T)/sin(w*l/a);u=0;forn=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;Page22雙曲型：波動(dòng)方程的求解與可視化figure(1)axis([01-55])h=plot(x,u(1,:),'linewidth',3);set(h,'erasemode','xor');forj=2:length(t);set(h,'ydata',u(j,:));axis([01-55])drawnowendfigure(2)waterfall(X(1:50:3000,:),T(1:50:3000,:),u(1:50:3000,:))Xlabel('x')Ylabel('t')Page23雙曲型：波動(dòng)方程的求解與可視化Page24圖2波動(dòng)方程解析解的分布例3長為l的桿，兩端自由，初始位移為c0x/l，初速度為零，研究其運(yùn)動(dòng)，定解問題為：Page25雙曲型：波動(dòng)方程的求解與可視化問題的解析解為：

其中：取C0=0.05，a=1，l=1，程序如下：Page26functiongzdN=50;t=0:0.005:2.0;x=0:0.002:1;ww=gzdfun(N,0);subplot(2,1,1)h1=plot(x,ww,'linewidth',3);set(h1,'erasemode','xor');axis([0,1,0,0.05]);xx=1:10:length(x);yy=0*xx;subplot(2,1,2)h2=plot(x(xx),yy,'r.','marker','.','markersize',25);set(h2,'erasemode','xor');axis([0,1.05,-0.1,0.1]);forn=2:length(t)ww=gzdfun(N,t(n));set(h1,'ydata',ww);uu=ww(xx)+x(xx);set(h2,'xdata',uu);drawnow;pause(0.02)endend雙曲型：波動(dòng)方程的求解與可視化functionwtx=gzdfun(N,t)x=0:0.002:1;a=1;wtx=1/2*0.05;fork=1:2:NBk=-4/(k*k*pi*pi)*cos(k*pi*t)*cos(k*pi*x)*0.05;wtx=wtx+Bk;endendPage27雙曲型：波動(dòng)方程的求解與可視化用y坐標(biāo)表示位移用質(zhì)點(diǎn)在水平方向上的移動(dòng)表示位移拋物型：熱傳導(dǎo)方程的求解與可視化熱傳導(dǎo)方程屬于拋物線方程，在MATLAB中是指如下形式：求解拋物線方程使用如下命令：u=parabolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)parabolic函數(shù)性質(zhì)與hyperbolic大致相同．Page28初始條件ut即是網(wǎng)格坐標(biāo)描述矩陣決定方程的類型時(shí)間矩陣邊界條件拋物型：熱傳導(dǎo)方程的求解與可視化例4求解下列熱傳導(dǎo)定解問題．求解域是方形區(qū)域，其中空間坐標(biāo)的個(gè)數(shù)由具體問題確定．Page29拋物型：熱傳導(dǎo)方程的求解與可視化【解】步驟如下：(1)題目定義g='squareg';%定義單位方形區(qū)域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.4);%ix是符合r<0.4的矩陣

u0(ix)=ones(size(ix));%產(chǎn)生行數(shù)與ix的行數(shù)相同的全1方陣(4)在時(shí)間段是0～0.1的20個(gè)點(diǎn)上求解

nframes=20;tlist=linspace(0,0.1,nframes);Page30拋物型：熱傳導(dǎo)方程的求解與可視化(5)求解此拋物問題

u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);(6)矩形網(wǎng)格插值x=linspace(-1,1,31);y=x;[unused,tn,a2,a3]=tri2grid(p,t,u0,x,y);(7)動(dòng)畫圖示結(jié)果

newplot;Mv=moviein(nframes);umax=max(max(u1));umin=min(min(u1));forj=1:nframes,...u=tri2grid(p,t,u1(:,j),tn,a2,a3);%用tri2grid的第三種形式，以最快速度插值Page31拋物型：熱傳導(dǎo)方程的求解與可視化i=find(isnan(u));%isnan(isnotanumber)當(dāng)不是數(shù)時(shí)返回1，是數(shù)時(shí)返回0，find則是找出非零元素

u(i)=zeros(size(i));

surf(x,y,u);%畫出由（x,y,z）決定的表面

caxis([uminumax]);colormap(cool),axis([-11-1101]);Mv(:,j)=getframe;endPage32拋物型：熱傳導(dǎo)方程的求解與可視化Page33圖3某一瞬時(shí)的熱傳導(dǎo)方程解分布拋物型：熱傳導(dǎo)方程的求解與可視化例

5討論如下的有限長細(xì)桿的熱傳導(dǎo)定解問題：Page34其中取l=40，a=20，且：拋物型：熱傳導(dǎo)方程的求解與可視化【解】定解問題的解析解（溫度分布）為：Page35

取將10個(gè)分波的合成，計(jì)算機(jī)仿真程序如下：functionyhjN=10;t=1e-5:0.00001:0.005;x=0:0.5:40;w=rcdf(N,t(1));h=plot(x,w,'linewidth',5);axis([0,40,0,1.5])forn=2:length(t)w=rcdf(N,t(n));set(h,'ydata',w);drawnow;end拋物型：熱傳導(dǎo)方程的求解與可視化Page36

functionu=rcdf(N,t)x=0:0.5:40;u=0;fork=1:2*Ncht=2/k/pi*(cos(k*pi*10/40.0)-cos(k*pi*30./40))*sin(k*pi*x./40);u=u+cht*exp(-(k^2*pi^2*20^2/1600*t));end圖4熱傳導(dǎo)細(xì)桿的溫度分布橢圓型：穩(wěn)定場方程的求解與可視化穩(wěn)定場方程屬于橢圓方程，下面我們來求解含有源項(xiàng)的標(biāo)準(zhǔn)穩(wěn)定場方程，也即是泊松方程．在MATLAB中是指如下形式：

求解橢圓方程用命令：u=assempde(b,p,e,t,c,a,f)

各輸入量的意義同前．由于是穩(wěn)定場方程，則輸出的矩陣u是坐標(biāo)矩陣p相應(yīng)的解．Page37橢圓型：穩(wěn)定場方程的求解與可視化例6用MATLAB求解下列泊松方程，并將計(jì)算機(jī)仿真解與精確解比較，方程如下：方法1：其解可以用偏微分方程工具箱求得，并繪制出其圖形。方法2：求出解析解，再利用MATLAB編程繪制出其解析解和誤差解的圖形分布。【解】滿足邊值條件的解析解（精確解）為：Page38橢圓型：穩(wěn)定場方程的求解與可視化計(jì)算機(jī)仿真：MATLAB求解步驟如下，并將仿真結(jié)果與精確解（解析解）進(jìn)行比較：(1)題目定義g='circleg';%單位圓

b='circleb1';%邊界零條件c=1;a=0;f=1; (2)初始的粗糙網(wǎng)格化

[p,e,t]=initmesh(g,'hmax',1);%hmax是內(nèi)部常數(shù)，每個(gè)三角形網(wǎng)格的大小不能超過hmax；1代表三角形網(wǎng)格個(gè)數(shù)增加的速度（一般在1.3左右）Page39橢圓型：穩(wěn)定場方程的求解與可視化(3)迭代直至得到誤差允許范圍內(nèi)的合理解error=[];er=1;whileer>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=[errorer];fprintf('Error:%e.Numberofnodes:%d\n',er,size(p,2));end%運(yùn)行結(jié)果是Error:1.292265e-002.Numberofnodes:25……Page40橢圓型：穩(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-exact,'zdata',u-exact,'mesh','on');%誤差解圖顯示Page41橢圓型：穩(wěn)定場方程的求解與可視化Page42圖5精確解與仿真解的誤差解圖(a)泊松方程的解析解圖(b)精確解與仿真解的誤差解橢圓型：穩(wěn)定場方程的求解與可視化

例7在矩形區(qū)域0<x<a，-b/2<y<b/2上，對滿足泊松方程

，且邊界上的值為零的定解問題的解，給出計(jì)算機(jī)仿真圖形。

【解】所討論的定解問題即為：Page43橢圓型：穩(wěn)定場方程的求解與可視化解析解的表達(dá)式：

計(jì)算機(jī)仿真程序：（程序中取a=8，b=8）symsaba=8;b=8;[X,Y]=meshgrid(0:0.2:a,-b/2:0.2:b/2)Z1=0;forn=1:1:10Z2=a^4*b*((-1)^n*n^2*pi^2+2-2*(-1)^n)*sinh(n*pi.*Y/a).*sin(n*pi.*X/a)/(n^5*pi^5*sinh(n*pi*b/(2*a)));Z1=Z1+Z2;end

Page44橢圓型：穩(wěn)定場方程的求解與可視化Z=Z1+X.*Y.*(a^3-X.^3)/12;colormap(hot);mesh(X,Y,Z)view(119,7);Page45圖6矩形域的泊松方程解的分布Page46本征值問題的求解與可視化二維本征值問題矩形區(qū)域的本征模與本征振動(dòng)邊長為b和c的的四周固定的矩形膜振動(dòng)的本征值問題為采用分離變量法可以得到本征模和本征值為Page47本征值問題的求解與可視化繪制前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.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);w31=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,w31);L=