课堂授课专题4：数学物理方程的计算机求解和可视化

数学物理建模与计算机辅助设计专题4：数学物理方程的计算机求解和可视化Page2本专题主要内容与参考资料

主要内容偏微分方程的计算机仿真求解方法双曲型（Hyperbolic）：波动方程的求解与可视化抛物型（Parabolic）：热传导方程的求解与可视化椭圆型（Elliptic）：稳定场方程的求解与可视化特征值问题的求解与可视化利用Matlab的PDE工具箱求解并进行可视化参考资料杨华军,数学物理方法,电子工业出版社彭芳麟,数学物理方程的Matlab解法与可视化,清华大学出版社仿真求解偏微分方程的类型通用的线性二阶偏微分方程：偏微分方程类型分为：双曲型方程：

抛物型方程：

椭圆型方程：

特征值问题：

特征值偏微分方程中不含参数f.Page3偏微分方程的仿真求解方法偏微分方程的计算机仿真求解方法：（1）MATLAB的偏微分方程工具箱（PDEToolbox）

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

典型偏微分的解的静态（或动态）三维可视化。Page4有限元法定义将连续的求解域离散成一组有限个、按照一定方式相互连结在一起的单元的组合体。将PDE转换成离散的线性代数方程系统进行求解。特点各种复杂单元可以用来对复杂的几何形状的求解域进行模型化处理。各节点上的解的近似函数可以用来求解整个求解域上任意点的结果。Page5MATLAB的偏微分方程工具箱（PDEToolbox）用图形用户界面（GraphicalUserInterface，简记作GUI）求解PDE问题主要采用以下三个步骤：设置PDE的定解问题．

即设置二维定解区域、边界条件以及方程的形式和系数；(2)用有限元法（FEM）求解PDE．即网格的生成、方程的离散以及求出数值解；

Mesh：生成网格，自动控制网格参数．Solve：求解．设置初始边值条件后，能给出t时刻的解；可以求出区间内的特征值．求解后可以加密网格再求解．(3)解的可视化．从GUI使用Plot方法实现可视化．用Color、Height、Vector等作图．对于含时方程，还可以生成解的动画．用PDEToolbox可以求解的基本方程有：椭圆方程、抛物方程、双曲方程、特征值方程、椭圆方程组以及非线性椭圆方程。．Page6偏微分方程工具箱求解定解问题例1解热传导方程

边界条件是齐次类型（

u=0），定解区域自定。【解】第一步：启动MATLAB，键入命令pdetool并回车，就进入GUI．

Options→Grid，打开栅格．第二步：选定定解区域绘制椭圆E1、圆E2、矩形R1、圆E3;在Setformula栏中键入E1-E2+R1-E3.第三步：选取边界Boundary→BoundaryMode，进入边界模式；Boundary→RemoveAllSubdomainBorders，去掉子域边界;Boundary→SpecifyBoundaryConditions

→BoundaryConditions→齐次Diriclet条件.Page7偏微分方程工具箱求解定解问题边界条件：解方程所需要的边界条件可以是以下两种形式：Page8狄利克里（Diriclet）边界条件广义诺依曼（Neumann）边界条件齐次边界条件：g=0，r=0第一类边界条件：Diriclet边界条件第三类边界条件：Neumann边界条件第二类边界条件：q=0偏微分方程工具箱求解定解问题第四步：设置方程类型PDEMode→PDESecification，选择方程类型:抛物型第五步：划分网格

Mesh→InitializeMesh，网格剖分；

Mesh→RefineMesh，网格密集化.第六步：解偏微分方程并显示图形解

Solve→SolvePDE，解偏微分方程并显示图形解。第七步：三维可视化

Plot

→Parameter，选择Color,Height(3-Dplot)，Showmesh第八步：绘制等值线图和矢量场图

Plot→Parameter，选择Contour和Arrows．PlotPage9双曲型：波动方程的求解与可视化双曲型：波动方程的求解与静态（或动态）三维可视化

1.求解双曲型方程求解双曲型方程调用形式如下：u1=hyperbolic（u0,ut0,tlist,b,p,e,t,c,a,f,d）Page10（a、c、d、f是参数）初始条件ut即是网格坐标描述矩阵决定方程的类型时间矩阵边界条件双曲型：波动方程的求解与可视化网格初始化命令：(1)[p,e,t]=initmesh(g)将求解区域进行三角形网格化，输出的p、e、t是网格数据．p-描述网格中点的x、y坐标．e-边缘矩阵，t-三角矩阵，描述区域的顶点．g-描述求解区域几何形状(2)[p,e,t]=refinemesh(g,p,e,t)

迭代过程，得到更细小的网格，使结果更精确．Page11双曲型：波动方程的求解与可视化2.动画图形显示

将所得的解形象地表示出来，为了加速绘图，首先把三角形网格转化成矩形网格．调用形式如下：

(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)用此命令之前，应先用一个tri2grid命令得出矩阵tn、a2、a3．用此方法可以加快速度．Page12三角形网格的矩阵矩形网格的坐标点各时刻三角形网格中的解插值法求得矩形网格点上的u值内插法的系数格点的指针矩阵双曲型：波动方程的求解与可视化主要的绘图（包括动画）命令函数有：moviein、movie、pedplot、pdesurf等Page13双曲型：波动方程的求解与可视化例1用MATLAB求解下面波动方程定解问题并动态显示解的分布

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

方法2：求出解析解，再利用MATLAB编程绘制出其解析的图形分布。Page14双曲型：波动方程的求解与可视化【解】采用步骤如下（1）题目定义

g='squareg';%定义单位方形区域

b='squareb3';%左右零边界条件，顶底零导数边界条件

c=1;a=0;f=0;d=1;（2）初始的粗糙网格化

[p,e,t]=initmesh('squareg');（3）初始条件

x=p(1,:)';%注意坐标向量都是列向量y=p(2,:)';u0=atan(sin(pi/2*x));ut0=2*cos(pi*x).*exp(cos(pi/2*y));Page15双曲型：波动方程的求解与可视化（4）在时间段0～5内的31个点上求解n=31;

tlist=linspace(0,5,n);%在0～5之间产生n个均匀的时间点（5）求解此双曲问题u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);得到如下结果：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双曲型：波动方程的求解与可视化解u1的动态图形显示：（6）矩形网格插值

delta=-1:0.1:1;[uxy,tn,a2,a3]=tri2grid(p,t,u1(:,1),delta,delta);gp=[tn;a2;a3];（7）在0～5时间内动画显示newplot;%建立新的坐标系M=moviein(n);umax=max(max(u1));umin=min(min(u1));Page17双曲型：波动方程的求解与可视化fori=1:nifrem(i,10)==0%当n是10的整数倍时，fprintf('%d’,i);%在命令窗口打印出相应的数字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双曲型：波动方程的求解与可视化%运行结果是：102030done动态解图可以直接通过MATLAB仿真程序执行看出，图1是动态图的某一瞬间的解的分布。Page19双曲型：波动方程的求解与可视化Page20图1某一瞬时的波动方程的解图双曲型：波动方程的求解与可视化Page21例2讨论弦的一端x=0,固定，x=L一端受迫作谐振动2sinωt，弦的初始位移和初始速度为零，给出弦振动的解图。【解】

根据题意得定解问题为：

解析解为：

双曲型：波动方程的求解与可视化计算机仿真程序如下：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双曲型：波动方程的求解与可视化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双曲型：波动方程的求解与可视化Page24图2波动方程解析解的分布例3长为l的杆，两端自由，初始位移为c0x/l，初速度为零，研究其运动，定解问题为：Page25双曲型：波动方程的求解与可视化问题的解析解为：

其中：取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双曲型：波动方程的求解与可视化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双曲型：波动方程的求解与可视化用y坐标表示位移用质点在水平方向上的移动表示位移抛物型：热传导方程的求解与可视化热传导方程属于抛物线方程，在MATLAB中是指如下形式：求解抛物线方程使用如下命令：u=parabolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)parabolic函数性质与hyperbolic大致相同．Page28初始条件ut即是网格坐标描述矩阵决定方程的类型时间矩阵边界条件抛物型：热传导方程的求解与可视化例4求解下列热传导定解问题．求解域是方形区域，其中空间坐标的个数由具体问题确定．Page29抛物型：热传导方程的求解与可视化【解】步骤如下：(1)题目定义g='squareg';%定义单位方形区域b='squareb1';%定义零边界条件c=1;a=0;f=1;d=1;(2)网格化[p,e,t]=initmesh(g);(3)定义初始条件

u0=zeros(size(p,2),1);%产生零矩阵u0，size(p,2)返回p的列数ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);%ix是符合r<0.4的矩阵

u0(ix)=ones(size(ix));%产生行数与ix的行数相同的全1方阵(4)在时间段是0～0.1的20个点上求解

nframes=20;tlist=linspace(0,0.1,nframes);Page30抛物型：热传导方程的求解与可视化(5)求解此抛物问题

u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);(6)矩形网格插值x=linspace(-1,1,31);y=x;[unused,tn,a2,a3]=tri2grid(p,t,u0,x,y);(7)动画图示结果

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抛物型：热传导方程的求解与可视化i=find(isnan(u));%isnan(isnotanumber)当不是数时返回1，是数时返回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抛物型：热传导方程的求解与可视化Page33图3某一瞬时的热传导方程解分布抛物型：热传导方程的求解与可视化例

5讨论如下的有限长细杆的热传导定解问题：Page34其中取l=40，a=20，且：抛物型：热传导方程的求解与可视化【解】定解问题的解析解（温度分布）为：Page35

取将10个分波的合成，计算机仿真程序如下：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抛物型：热传导方程的求解与可视化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热传导细杆的温度分布椭圆型：稳定场方程的求解与可视化稳定场方程属于椭圆方程，下面我们来求解含有源项的标准稳定场方程，也即是泊松方程．在MATLAB中是指如下形式：

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

各输入量的意义同前．由于是稳定场方程，则输出的矩阵u是坐标矩阵p相应的解．Page37椭圆型：稳定场方程的求解与可视化例6用MATLAB求解下列泊松方程，并将计算机仿真解与精确解比较，方程如下：方法1：其解可以用偏微分方程工具箱求得，并绘制出其图形。方法2：求出解析解，再利用MATLAB编程绘制出其解析解和误差解的图形分布。【解】满足边值条件的解析解（精确解）为：Page38椭圆型：稳定场方程的求解与可视化计算机仿真：MATLAB求解步骤如下，并将仿真结果与精确解（解析解）进行比较：(1)题目定义g='circleg';%单位圆

b='circleb1';%边界零条件c=1;a=0;f=1; (2)初始的粗糙网格化

[p,e,t]=initmesh(g,'hmax',1);%hmax是内部常数，每个三角形网格的大小不能超过hmax；1代表三角形网格个数增加的速度（一般在1.3左右）Page39椭圆型：稳定场方程的求解与可视化(3)迭代直至得到误差允许范围内的合理解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%运行结果是Error:1.292265e-002.Numberofnodes:25……Page40椭圆型：稳定场方程的求解与可视化(4)把结果用图形表示：%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椭圆型：稳定场方程的求解与可视化Page42图5精确解与仿真解的误差解图(a)泊松方程的解析解图(b)精确解与仿真解的误差解椭圆型：稳定场方程的求解与可视化

例7在矩形区域0<x<a，-b/2<y<b/2上，对满足泊松方程

，且边界上的值为零的定解问题的解，给出计算机仿真图形。

【解】所讨论的定解问题即为：Page43椭圆型：稳定场方程的求解与可视化解析解的表达式：

计算机仿真程序：（程序中取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椭圆型：稳定场方程的求解与可视化Z=Z1+X.*Y.*(a^3-X.^3)/12;colormap(hot);mesh(X,Y,Z)view(119,7);Page45图6矩形域的泊松方程解的分布Page46本征值问题的求解与可视化二维本征值问题矩形区域的本征模与本征振动边长为b和c的的四周固定的矩形膜振动的本征值问题为采用分离变量法可以得到本征模和本征值为Page47本征值问题的求解与可视化绘制前4个本征函数的图形%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个本征函数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=