matlab描述几何区域的文件,应用Matlab的PDE Toolbox求解偏微分方程综合实验

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

2011 届信计专业学生综合实验题目(要求按照所附开题报告表的格式填写提交开题报告。一个小组选做一题，小组全体成员共同完成，每个小组只提交一份实验报告，按照出力多少排名。提交时间在本学期 18周以前。 )3 应用 Matlab 的 PDE Toolbox 求解偏微分方程熟悉 Matlab 的 PDE 工具箱的功能，并用其求解具有工程背景的偏微分方程，要求分别对三种类型方程：抛物型、椭圆形和双曲型。阐述清楚如下方法：1、这里我们先脱离问题所含有的工程背景，分别举三个例子，大致的说明一下如何使用 Matlab的 PDE 工具箱来求解三种类型的偏微分方程。(1) 椭圆形方程： infauc,)(考虑一块圆形金属片，中心挖去一正方形，外边界满足 Neumann 条件，内边界满足Dirichlet 条件： .),(xikeyvr考虑到入射波以 方向，所以上式可以写成 这样得到求解这个入射波x,ikxer的定解问题：),(yxr外 边 界 上内 边 界 上在 区 域 内,60,2iiknrexiix这里取 波长 为 0.1。,60k现在用 GUI 来求解上述这个问题，并最终获得其解得图形：图 1 初始网格、加密以及网格剖分数据图 2 解的三维图形图 3 解的二维动画图附录 1：二维动画的 Matlab 程序echo onclc%程序段一：解 Helmholtz 方程-div(grad(u))-k^2u=0%并研究正方形上的反射波，波源来自右边clc%入射波波数为 60k=60;g= scatterg ;%scatterg 描述几何区域的文件名，此区域为园内有一方洞b= scatterb ;%scatterb 是描述边界条件的文件名，内边界满足%Dirichlet 条件，外边界满足 Neumann 条件%选择方程系数 c,a,fc=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 equalclc%在复平面上求解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 offM=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(2) 抛物型方程： infauctud,)(考虑一个圆柱形放射性杆，其左端供热，右端保持常温，侧面与环境有热交换。由于放性作用，热量均匀地产生。初始温度为 。于是可以用如下方程描述： C0fuktc)(其中 为密度， 为杆的热容量， 为导热系数， 为放射性热源密度。这一金属杆的密度c取为 热容量 为 导热系数 为 热源密度,/7803mkgc,/500Ckgwsk,/40Cmw为 右端恒温为 侧面环境温度为 热交换系数为f2w,1,1左端的热流为,/50C./2边界条件(如右图 4)：在杆的左端( ):处5.1zrucn50)(在杆的右端( ):处在杆的侧面( )： 处2.0r 10)(r在杆的轴心( ):处0ucn初始温度： |0tu现在用 GUI 来求解上述这个问题，并最终获得其解得图形：图 5 解的 2 维动画图形 图 6 解的 3 维动画图形图 7 解的动画图形比较图(3) 双曲型方程： infauctud,)(2考虑如下二维波动方程的定界问题，并最终获得其解得图形：  )2sin(12 i3),2(arctnos,0|0| 1,1|),(),( yyx extuxutu yxyyxt 现在用 GUI 来求解上述这个问题，并最终获得其解得图形：图 8 解的 3 维动画图形附录 2：双曲线三维动画的 Matlab 程序[p.e.t]=initmesh( squareg );%初始化网格x=p(1,:) ,y=p(2,:) ;u0=atan(cos(pi/2*x));%ut0=3*sin(pi*x).*exp(sin(pi/2*y));%n=31;tlist=linspace(0,5,n);%uu=hyperbolic(u0,ut0,tlist,’squareb3’,p,e,t,1,0,0,1);%delta=-1:0.1:1;[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),detla,delta);%gp=[tn;a2;a3];umax=max(max(uu));umin=min(min(uu));newplotM=moviein(n);for i=1:n,pdeplot(p,e,t,’xydata’,uu(:,i),’zdata’,uu(:,i),…‘mesh’,’off’,’xygrid’,’on’,’gridparam’,gp,…‘colorbar’,’off’,’zstyle’,’continuous’);axis([-1 1 -1 1 umin umax]);caxis([umin umax]);M(:,i)=getframe;endMovie(M,10);2、向问题中添加一定的工程背景，以下我们详尽的说明如何使用 Matlab 的 PDE 工具箱来求解三种类型的偏微分方程。(1)椭圆形方程： infauc,)(在磁铁、电动机、变压器这类问题中，考虑如下定常的 Maxwell 方程：HB0,J其中 B 是磁感应强度，H 是磁场强度，J 是电流密度， 是材料的磁导率。由于故存在一个静磁矢势 A，使得0,及.)1(JA平面问题中假设电流平行于 z 轴，故 A 仅有 z 分量 ),0(),0(J从而，上述方程可以简化为椭圆型 PDE：J