电磁场与电磁波研究性学习报告—有限差分方法

1、研究型学习报告 利用有限差分法解决静电场第一类边值问题马啸13274015理论分析理论分析22220|( )Lxyf sxyoL10234hhD D理论分析理论分析 有限差分法拉普拉斯方程的离散化 将区域D划分为正方形网格，在节点处代入导数 公式，可得 将其带入拉普拉斯方程，可得方程的差分格式为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 迭代法 分量形式 向量形式 其中 参考文献给出，最佳松弛因子的计算公式为 1111111,1,2,., ;0,1,2,.inkkkkiiiijjijjjj iiixxba xa xin ka 1kkxB xf11BDLDU1fDLb2opt1+ sin h理论分

3、析理论分析计算机仿真计算机仿真 问题描述 如图所示，有一长方形的导体槽，a = 20，b = 5，设槽的长度为无限长，槽上有两块与槽绝缘的盖板，电位分别为100V和 50V，其它板电位为0V，求槽内的电位分布.yx10050ab计算机仿真计算机仿真 获取系数矩阵%get_coeeficient_matrix_function%输入未知数矩阵的行数输入未知数矩阵的行数m及列数及列数n,根据有限差分法返根据有限差分法返回方程的系数矩阵及常数列向量回方程的系数矩阵及常数列向量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;%将矩阵个数的未知数转换为自然数排将矩阵个数的未知数转换为自然数排列列 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;%第一行完毕第一行完毕计算机仿真计算机仿真 获取系数

5、矩阵B=0.25.*B;%得到得到Jacobi迭代矩阵迭代矩阵 for i=1:m*n B(i,i)=-1;endA=-B;%得到系数矩阵得到系数矩阵 b=b/4;%A=(A-D)+DD=zeros(m*n,m*n);%A所对应的对角阵所对应的对角阵for i=1:m*n D(i,i)=A(i,i);endf=D*b;b=f;计算机仿真计算机仿真 Gauss-Seidel迭代法h=input(所期望步长所期望步长:);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); %创建符合步长要求的

6、创建符合步长要求的矩阵存放各标志点电位值矩阵存放各标志点电位值 for j=2:col_h-1 for i=2:line_h-1 u(i,j)=50+50/line_h*(i-1); %初始初始条件按内部线性变化确定条件按内部线性变化确定end %初始条件对结果影响不大，但影响循环次数初始条件对结果影响不大，但影响循环次数end %完成差分方程初始条件的设定完成差分方程初始条件的设定A,b=get_coefficient_matrix_function(line_x,col_x);%使使用自定义函数获得方程组系数矩阵及常数列向量用自定义函数获得方程组系数矩阵及常数列向量mat,flag=cho

7、l(A);%检查检查A是否正正定是否正正定if flag=0 disp(收敛收敛)else disp(不收敛不收敛)end %若若A正定，可以使用正定，可以使用Gauss-Seidel迭代法求解方程迭代法求解方程计算机仿真计算机仿真while flag=0 %在迭代法收敛的条件下进入循环在迭代法收敛的条件下进入循环 while (t1e-6) n=n+1; t=0; %必须将误差控制项置零必须将误差控制项置零 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);%用无穷大范数作为误差控制用无穷大范数作为误差控制 end end end u=v; endend计算机仿真计算机仿真 Successive Over Relaxation 迭代法w_sor=2/(1+sin(pi*h);%计算最佳松弛因子计算最佳松弛因子w_sor;%打印松弛因子以备查验打印松弛因子以备查验if (1w_sor)&(w_sor2) %在迭代法收敛的条件下进入循环在迭代法收敛的条件下进入循环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迭代法迭代法分量形式分量形式计算机仿真计算机仿真 仿真结果计算机仿真计算机仿真 程序代码for i=2:line_h-1 for j=2:col_h-1 x=(j-1)*h; y=(i-1)*h;%将矩阵位置转化为直角坐将矩阵位置转化为直角坐标标 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为解析解表达式，手工算出为解析解表达式，手工算出 end u_analysis(i,j)=phi; error(i