应用Matlab的PDE Toolbox求解偏微分方程综合实验.doc

PAGE 2011届信计专业学生综合实验题目 （要求按照所附开题报告表的格式填写提交开题报告。一个小组选做一题，小组全体成员共同完成，每个小组只提交一份实验报告，按照出力多少排名。提交时间在本学期18周以前。） 3 应用Matlab的PDE Toolbox求解偏微分方程 熟悉Matlab的PDE工具箱的功能，并用其求解具有工程背景的偏微分方程，要求分别对三种类型方程：抛物型、椭圆形和双曲型。阐述清楚如下方法： 1、这里我们先脱离问题所含有的工程背景，分别举三个例子，大致的说明一下如何使用Matlab的PDE工具箱来求解三种类型的偏微分方程。 椭圆形方程： 考虑一块圆形金属片，中心挖去一正方形，外边界满足Neumann条件，内边界满足Dirichlet条件： 考虑到入射波以方向，所以上式可以写成这样得到求解这个入射波的定解问题： 这里取波长为0.1。 现在用GUI来求解上述这个问题，并最终获得其解得图形： 图1 初始网格、加密以及网格剖分数据 图2 解的三维图形 图3 解的二维动画图 附录1：二维动画的Matlab程序 echo on clc %程序段一：解Helmholtz方程-div(grad(u))-k^2u=0 %并研究正方形上的反射波，波源来自右边 clc %入射波波数为60 k=60; g=scatterg;%scatterg描述几何区域的文件名，此区域为园内有一方洞 b=scatterb;%scatterb是描述边界条件的文件名，内边界满足 %Dirichlet条件，外边界满足Neumann条件 %选择方程系数c,a,f c=1; a=-k^2; f=0; %程序段二：初始化网格和加密网格 [p,e,t]=initmesh(g); [p,e,t]=refinemesh(g,p,e,t); [p,e,t]=refinemesh(g,p,e,t); %会出网格图 pdemesh(p,e,t);axis equal clc %在复平面上求解 u=assempde(b,p,e,t,c,a,f); %取复数解的实部 h=newplot;set(get(h,Parent),Renderer,zbuffer) pdeplot(p,e,t,xydata,real(u),zdata,real(u),... mesh,off); colormap(cool) clc %制作反射波的动画程序 m=10;%帧数 h=newplot;hf=get(h,Parent);set(hf,Renderer,zbuffer) axis tight,set(gca,DataAspectRatio,[1 1 1]); axis off M=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(done\n);... end,... end %显示动画 movie(hf,M,50); echo off 抛物型方程： 考虑一个圆柱形放射性杆，其左端供热，右端保持常温，侧面与环境有热交换。由于放 性作用，热量均匀地产生。初始温度为。于是可以用如下方程描述： 其中为密度，为杆的热容量，为导热系数，为放射性热源密度。这一金属杆的密度 取为热容量为导热系数为热源密度为右端恒温为侧面环境温度为热交换系数为左端的热流为 边界条件(如右图4)： 在杆的左端(): 在杆的右端(): 在杆的侧面()： 在杆的轴心(): 初始温度： 现在用GUI来求解上述这个问题，并最终获得其解得图形： 图5 解的2维动画图形 图6 解的3维动画图形 图7 解的动画图形比较图 双曲型方程： 考虑如下二维波动方程的定界问题，并最终获得其解得图形： 现在用GUI来求解上述这个问题，并最终获得其解得图形： 图8 解的3维动