有限元地MATLAB解法

1、标准 MATLAB1.打开。2.输入“”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是 partial differential equations的话可点击Options菜单下Grid命令，打开栅格。3.完成平面几何模型：在 PDE Toolbox的窗口中，点击工具栏下的 圆 Setformula栏进行编辑并（如双脊波导 R1+R2+R3改为 RI-R2-R3,设定、s/a、d/b的值从而方便下步设定坐标）用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。4.用左键双击矩形进行坐标设置：将大的矩形left和 bottom都设为 ，width是矩形波导

2、的 X轴的长度，height是矩形波导的 y轴5.进行边界设置：点击“Boundary”中的“Boundary ”,再点文案标准击“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。6.进入PDEPDEPDEMode”命令，进入PDE模PDESpecificationElliptic”为椭ParabolicHyperbolicEigenmodes”为特征值问题。7.InitializeMesh”进行初次剖分，若要剖的更细，再点击“Refine ”进行网格加密。8.

3、SolveSolvePDE解偏微分方程并显示图形解，u值即为Hz或者Ez。9.ParametersPlotSelection”对话框。选中Color,Height(3-D plot)和 Show mesh三项，然后单击“Plot”按钮，显示三维图形解。10.PlotParameters”选项打开“Plot ”对话框。选中 Contour和 Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。11.将计算结果条件和边界导入MATLABExportSolution,再点击“”中“Export 文案标准12.在 MATLAB中将编好的计算程序导入，按F5运行。备注：Property

4、（属性）用于画图时选用相应的绘图类型u方程的解abs(grad(u)abs(c*grad(u)- grad(u)每个三角形的中心的u的绝对值每个三角形的中心的u的绝对值u的负梯度u我们也可以用MATLAB程序求解PDE问题，同时显示解的图形；文案标准位为，求槽内的电位分布：0V0V（）标出各剖分点坐标值；（）画出等电位图。）编写以下程序得：x=0:5y=0:5文案标准X,Y=meshgrid(x,y)plot(X,Y)hold onplot(Y,X)for i=0:5s=i:5t=0:(5-i)plot(s,t)plot(t,s)end得到剖分图如下：（）用有限元法编写程序如下：文案标准Nx=

5、6;Ny=6;Xm=5;Ym=15;Np=5;Nq=5;for i=1:Nxfor j=1:NyN(i,j)=(i-1)*Ny+j; /i列行的节点编号/X(N(i,j)=(i-1)*Xm/Np;/节点横坐标/Y(N(i,j)=(j-1)*Ym/Nq;/节点纵坐标/endendfor i=1:2*Xmfor j=1:Ymif rem(i,2)=1L(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)=0L(i,j)=(i-1)*Ny+j;p(i,j)=(2i-2)*N

6、y/2+j;q(i,j)=p(i,j)+Ny;文案标准r(i,j)=q(i,j)+1;endendendfor i=1:2*Xmfor 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)=X(r(i,j)-X(q(i,j);c(q(i,j)=X(p(i,j)-X(r(i,j);c(r(i,j)=X(q(i,j)-X(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(Nx*Ny);

7、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(

8、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*Xmfor 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

9、,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:11K(i,:)=0;K(i,i)=1;endfor i=1:11:111K(i,:)=0;K(i,i)=1;endfor i=111:121K(i,:)=0;K(i,i)=1;end文案标准for i=11:11:121K(i,:)=

10、0;K(i,i)=1;endB=zeros(121,1);for i=11:11:121B(i,1)=100;endU=KB;b=1;XX=zeros(11,11)for j=1:11for i=1:11XX(i,j)=U(b,1);b=b+1;endendsubplot(1,2,1),mesh(XX)axis(0,11,0,11,0,100)subplot(1,2,2),contour(XX,15)文案标准hold on（）由上面的程序得到节点电位：V1=000000 7.1429 9.8214 7.1429 00 18.7500 25.0000 18.7500 00 42.8571 52.

11、6786 42.8571 00 100.0000 100.0000 100.0000 0（）由程序得到的电场分布图及等位线图如下：文案标准4.用有限元法求矩形波导（b/a=0.45）的：（）电场分布图；（2）求 TE 模式下的主模、第一、二高次模的截止波长（5（3）求 TM 模式下的主模、第一、二高次模的截止波长（5 MATLAB 中的 PDE ，窄边尺寸为 0.45a。（） 主模的电场分布图如下：在Neumann 边界条件下：文案标准在 Dirichlet边界条件下：（）在 TE 模式下设置边界条件为 Neumann 条件，使用编制好的程序计算出主模的截止波长为1.9988a，第一高次模为0.9977a，第二高次模为 0.8972a，截止波长图如下：a00.9977a时1.9988a时0.9977a1.9988a时 区10（）在 TM 模式下设置边界条件为 Dirichlet 条件，使用编制好的程序计算出主模的截止波长为 0.8179a，第一高次模为文案标准0.6655a，第二高次模为 0.5315a，截止波长图如下：3113a0时时时 区115.用时域有限差分求解上述 4题中的前两问。1）根据时域有限差分编写的程