有限元地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 ”進行網(wǎng)格加密。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列行的節(jié)點編號/X(N(i,j)=(i-1)*Xm/Np;/節(jié)點橫坐標/Y(N(i,j)=(j-1)*Ym/Nq;/節(jié)點縱坐標/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（）由上面的程序得到節(jié)點電位：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時 區(qū)10（）在 TM 模式下設置邊界條件為 Dirichlet 條件，使用編制好的程序計算出主模的截止波長為 0.8179a，第一高次模為文案標準0.6655a，第二高次模為 0.5315a，截止波長圖如下：3113a0時時時 區(qū)115.用時域有限差分求解上述 4題中的前兩問。1）根據(jù)時域有限差分編寫的程