電磁場(chǎng)與電磁波研究性學(xué)習(xí)報(bào)告—有限差分方法

1、研究型學(xué)習(xí)報(bào)告 利用有限差分法解決靜電場(chǎng)第一類邊值問(wèn)題馬嘯13274015理論分析理論分析22220|( )Lxyf sxyoL10234hhD D理論分析理論分析 有限差分法拉普拉斯方程的離散化 將區(qū)域D劃分為正方形網(wǎng)格，在節(jié)點(diǎn)處代入導(dǎo)數(shù) 公式，可得 將其帶入拉普拉斯方程，可得方程的差分格式為0021032222042222x xy yxhyhxyoL10234hhD D0404321理論分析理論分析 差分方程組的解Jacobi迭代法 分量格式 向量格式 其中 111,1,2,., ;0,1,2,.nkkiiijjjiij kxba xin kaJJxB xf11,JJBID A fD b

2、理論分析理論分析 差分方程組的解Gauss-Seidel 迭代法 分量形式 向量形式 其中 111111,1,2,., ;0,1,2,.inkkkiiijjijjjj iiixba xa xin ka 1kkSSxB xf1SB = DLU1SfDLb理論分析理論分析 差分方程組的解Successive Over Relaxation 迭代法 分量形式 向量形式 其中 參考文獻(xiàn)給出，最佳松弛因子的計(jì)算公式為 1111111,1,2,., ;0,1,2,.inkkkkiiiijjijjjj iiixxba xa xin ka 1kkxB xf11BDLDU1fDLb2opt1+ sin h理論分

3、析理論分析計(jì)算機(jī)仿真計(jì)算機(jī)仿真 問(wèn)題描述 如圖所示，有一長(zhǎng)方形的導(dǎo)體槽，a = 20，b = 5，設(shè)槽的長(zhǎng)度為無(wú)限長(zhǎng)，槽上有兩塊與槽絕緣的蓋板，電位分別為100V和 50V，其它板電位為0V，求槽內(nèi)的電位分布.yx10050ab計(jì)算機(jī)仿真計(jì)算機(jī)仿真 獲取系數(shù)矩陣%get_coeeficient_matrix_function%輸入未知數(shù)矩陣的行數(shù)輸入未知數(shù)矩陣的行數(shù)m及列數(shù)及列數(shù)n,根據(jù)有限差分法返根據(jù)有限差分法返回方程的系數(shù)矩陣及常數(shù)列向量回方程的系數(shù)矩陣及常數(shù)列向量function A,b=get_coefficient_matrix_function(m,n)B=zeros(m*n,m*

4、n);b=zeros(m*n,1);rup=50;rlw=100;rlf=0;rrt=0;i=1; j=1; k=(i-1)*n+j;%將矩陣個(gè)數(shù)的未知數(shù)轉(zhuǎn)換為自然數(shù)排將矩陣個(gè)數(shù)的未知數(shù)轉(zhuǎn)換為自然數(shù)排列列 B(k,k+1)=1; B(k,k+n)=1; b(k,1)=rup+rlf; for j=2:n-1 k=(i-1)*n+j; B(k,k-1)=1;B(k,k+1)=1;B(k,k+n)=1; b(k,1)=rup; end j=n; k=(i-1)*n+j; B(k,k-1)=1;B(k,k+n)=1; b(k,1)=rup+rlf;%第一行完畢第一行完畢計(jì)算機(jī)仿真計(jì)算機(jī)仿真 獲取系數(shù)

5、矩陣B=0.25.*B;%得到得到Jacobi迭代矩陣迭代矩陣 for i=1:m*n B(i,i)=-1;endA=-B;%得到系數(shù)矩陣得到系數(shù)矩陣 b=b/4;%A=(A-D)+DD=zeros(m*n,m*n);%A所對(duì)應(yīng)的對(duì)角陣所對(duì)應(yīng)的對(duì)角陣for i=1:m*n D(i,i)=A(i,i);endf=D*b;b=f;計(jì)算機(jī)仿真計(jì)算機(jī)仿真 Gauss-Seidel迭代法h=input(所期望步長(zhǎng)所期望步長(zhǎng):);col_h=20/h+1;col_x=col_h-2;line_h=5/h+1;line_x=line_h-2;u=zeros(line_h,col_h); %創(chuàng)建符合步長(zhǎng)要求的

6、創(chuàng)建符合步長(zhǎng)要求的矩陣存放各標(biāo)志點(diǎn)電位值矩陣存放各標(biāo)志點(diǎn)電位值 for j=2:col_h-1 for i=2:line_h-1 u(i,j)=50+50/line_h*(i-1); %初始初始條件按內(nèi)部線性變化確定條件按內(nèi)部線性變化確定end %初始條件對(duì)結(jié)果影響不大，但影響循環(huán)次數(shù)初始條件對(duì)結(jié)果影響不大，但影響循環(huán)次數(shù)end %完成差分方程初始條件的設(shè)定完成差分方程初始條件的設(shè)定A,b=get_coefficient_matrix_function(line_x,col_x);%使使用自定義函數(shù)獲得方程組系數(shù)矩陣及常數(shù)列向量用自定義函數(shù)獲得方程組系數(shù)矩陣及常數(shù)列向量mat,flag=cho

7、l(A);%檢查檢查A是否正正定是否正正定if flag=0 disp(收斂收斂)else disp(不收斂不收斂)end %若若A正定，可以使用正定，可以使用Gauss-Seidel迭代法求解方程迭代法求解方程計(jì)算機(jī)仿真計(jì)算機(jī)仿真while flag=0 %在迭代法收斂的條件下進(jìn)入循環(huán)在迭代法收斂的條件下進(jìn)入循環(huán) while (t1e-6) n=n+1; t=0; %必須將誤差控制項(xiàng)置零必須將誤差控制項(xiàng)置零 for i=2:line_h-1; for j=2:col_h-1; v(i,j)=(u(i,j+1)+u(i+1,j)+v(i-1,j)+v(i,j-1)/4; %Seidel迭代法迭

8、代法 if tabs(v(i,j)-u(i,j) t=abs(v(i,j)-u(i,j);%用無(wú)窮大范數(shù)作為誤差控制用無(wú)窮大范數(shù)作為誤差控制 end end end u=v; endend計(jì)算機(jī)仿真計(jì)算機(jī)仿真 Successive Over Relaxation 迭代法w_sor=2/(1+sin(pi*h);%計(jì)算最佳松弛因子計(jì)算最佳松弛因子w_sor;%打印松弛因子以備查驗(yàn)打印松弛因子以備查驗(yàn)if (1w_sor)&(w_sor2) %在迭代法收斂的條件下進(jìn)入循環(huán)在迭代法收斂的條件下進(jìn)入循環(huán)v(i,j)=u(i,j)+(u(i,j+1)+u(i+1,j)+v(i-1,j)+v(i,

9、j-1)-4*u(i,j)*w_sor/4; %SOR迭代法迭代法分量形式分量形式計(jì)算機(jī)仿真計(jì)算機(jī)仿真 仿真結(jié)果計(jì)算機(jī)仿真計(jì)算機(jī)仿真 程序代碼for i=2:line_h-1 for j=2:col_h-1 x=(j-1)*h; y=(i-1)*h;%將矩陣位置轉(zhuǎn)化為直角坐將矩陣位置轉(zhuǎn)化為直角坐標(biāo)標(biāo) phi=0; for m=1:2:M Cm=200/(m*pi); Bm=(400/(m*pi)-Cm*cosh(m*pi/4)/sinh(m*pi/4); phi_i=(Bm*sinh(m*pi*y/20)+Cm*cosh(m*pi*y/20)*sin(m*pi*x/20); phi=phi+phi_i; %phi為解析解表達(dá)式，手工算出為解析解表達(dá)式，手工算出 end u_analysis(i,j)=phi; error(i