有限元的MATLAB解法.doc
-有限元的MATLAB解法1.翻开MATLAB。2.输入“pdetool再回车，会跳出PDE Toolbo*的窗口PDE意为偏微分方程，是partial differential equations的缩写，需要的话可点击Options菜单下Grid命令，翻开栅格。3.完成平面几何模型：在PDE Toolbo*的窗口中，点击工具栏下的矩形几何模型进展制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进展编辑并如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标用算术运算符将图形对象名称连接起来，假设还需要，可进展储存，形成M文件。4.用左键双击矩形进展坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的*轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。5.进展边界设置：点击“Boundary中的“Boundary Mode,再点击“Boundary中的“Specify Boundary Conditions,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。6.进入PDE模式：点击"PDE"菜单下“PDE Mode命令，进入PDE模式，单击“PDE Specification,设置方程类型，“Elliptic为椭圆型，“Parabolic为抛物型，“Hyperbolic为双曲型，“Eigenmodes为特征值问题。7.对模型进展剖分：点击“Mesh中“Initialize Mesh进展初次剖分，假设要剖的更细，再点击“Refine Mesh进展网格加密。8.进展计算：点击“Solve中“Solve PDE,解偏微分方程并显示图形解，u值即为Hz或者Ez。9.单击“Plot菜单下“Parameters选项，翻开“Plot Selection对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单击“Plot按钮，显示三维图形解。10.如果要画等值线图和矢量场图，单击“Plot菜单下“Parameters选项,翻开“Plot Selection对话框。选中Contour和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。11.将计算结果条件和边界导入MATLAB中：点击“E*port Solution,再点击“Mesh中“E*port Mesh。12.在MATLAB中将编好的计算程序导入，按F5运行。备注：Property属性用于画图时选用相应的绘图类型u 方程的解abs(grad(u) 每个三角形的中心的u的绝对值abs(c*grad(u) 每个三角形的中心的c·u的绝对值- grad(u) u的负梯度-u我们也可以用MATLAB程序求解PDE问题，同时显示解的图形；一个长直接接地金属矩形槽，其侧壁与底面电位均为0，顶盖电位为100V，求槽的电位分布：1画出剖分图尺寸与书上一样；2标出各剖分点坐标值；3求出各点电位值用有限差分；4画出等电位图。解：1编写以下程序得：*=0:5y=0:5*,Y=meshgrid(*,y)plot(*,Y)hold onplot(Y,*)for i=0:5 s=i:5 t=0:(5-i) plot(s,t) plot(t,s)end得到剖分图如下：2用有限元法编写程序如下：N*=6;Ny=6;*m=5;Ym=15;Np=5;Nq=5;for i=1:N*for j=1:Ny N(i,j)=(i-1)*Ny+j; /i列j行的节点编号/ *(N(i,j)=(i-1)*m/Np;/节点横坐标/ Y(N(i,j)=(j-1)*Ym/Nq;/节点纵坐标/endendfor i=1:2*m for j=1:Ym if rem(i,2)=1 L(i,j)=(i-1)*Nq+j; p(i,j)=2*(i-1)*Ny/2+Ny+j+1; q(i,j)=p(i,j)-Ny; r(i,j)=q(i,j)-1; else rem(i,2)=0 L(i,j)=(i-1)*Ny+j; p(i,j)=(2i-2)*Ny/2+j; q(i,j)=p(i,j)+Ny; r(i,j)=q(i,j)+1; end endendfor i=1:2*m for j=1:Ymb(p(i,j)=Y(q(i,j)-Y(r(i,j);b(q(i,j)=Y(r(i,j)-Y(p(i,j);b(r(i,j)=Y(p(i,j)-Y(q(i,j);c(p(i,j)=*(r(i,j)-*(q(i,j);c(q(i,j)=*(p(i,j)-*(r(i,j);c(r(i,j)=*(q(i,j)-*(p(i,j);area(i,j)=(b(p(i,j)*c(q(i,j)-b(q(i,j)*c(p(i,j)/2; K=zeros(N*Ny); Kpp(i,j)=(b(p(i,j)2+c(p(i,j)2)/(2*area(i,j);Kpq(i,j)=(b(p(i,j)*b(q(i,j)+c(p(i,j)*c(q(i,j)/(2*area(i,j); Kpr(i,j)=(b(p(i,j)*b(r(i,j)+c(p(i,j)*c(r(i,j)/(2*area(i,j);Kqp(i,j)=Kpq(i,j); Kqq(i,j)=(b(q(i,j)2+c(q(i,j)2)/(2*area(i,j);Kqr(i,j)=(b(q(i,j)*b(r(i,j)+c(q(i,j)*c(r(i,j)/(2*area(i,j);Krp(i,j)=Kpr(i,j); Krq(i,j)=Kqr(i,j); Krr(i,j)=(b(r(i,j)2+c(r(i,j)2)/(2*area(i,j); endendfor i=1:2*m for j=1:YmK(p(i,j),p(i,j)=Kpp(i,j)+K(p(i,j),p(i,j);K(p(i,j),q(i,j)=Kpq(i,j)+K(p(i,j),q(i,j);K(p(i,j),r(i,j)=Kpr(i,j)+K(p(i,j),r(i,j);K(q(i,j),p(i,j)=Kqp(i,j)+K(q(i,j),p(i,j);K(q(i,j),q(i,j)=Kqq(i,j)+K(q(i,j),q(i,j);K(q(i,j),r(i,j)=Kqr(i,j)+K(q(i,j),r(i,j);K(r(i,j),p(i,j)=Krp(i,j)+K(r(i,j),p(i,j);K(r(i,j),q(i,j)=Krq(i,j)+K(r(i,j),q(i,j);K(r(i,j),r(i,j)=Krr(i,j)+K(r(i,j),r(i,j); endendfor i=1:11 K(i,:)=0; K(i,i)=1;endfor i=1:11:111 K(i,:)=0; K(i,i)=1;endfor i=111:121 K(i,:)=0; K(i,i)=1;endfor i=11:11:121 K(i,:)=0; K(i,i)=1;endB=zeros(121,1);for i=11:11:121 B(i,1)=100;endU=KB;b=1;*=zeros(11,11)for j=1:11 for i=1:11 *(i,j)=U(b,1); b=b+1; endendsubplot(1,2,1),mesh(*)a*is(0,11,0,11,0,100)subplot(1,2,2),contour(*,15)hold on3由上面的程序得到节点电位：4由程序得到的电场分布图及等位线图如下：. z.-. z.-4.用有限元法求矩形波导b/a=0.45的： 1电场分布图； 2求TE模式下的主模、第一、二高次模的截止波长5次，画出截至波长图； 3求TM模式下的主模、第一、二高次模的截止波长5次，画出截至波长图。解：利用MATLAB中的PDE工具箱：取矩形波导的宽边尺寸为a，窄边尺寸为0.45a。（1） 主模的电场分布图如下： 在Neumann边界条件下：在Dirichlet边界条件下：2在TE模式下设置边界条件为Neumann条件，使用编制好的程序计算出主模的截止波长为1.9988a，第一高次模为0.9977a，第二高次模为0.8972a，截止波长图如下： 3在TM模式下设置边界条件为Dirichlet条件，使用编制好的程序计算出主模的截止波长为0.8179a，第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下：5.用时域有限差分求解上述4题中的前两问。解：1根据时域有限差分编写的程序可画出主模的电场分布图如下： 在Neumann边界条件下：在Dirichlet边界条件下：2根据时域有限差分编写的程序可画出频谱图和场构造图，从左图中可以读出主模截止频率值，主模，根据，其中, 从而计算出主模截至波长。. z.
