zuoye.m

来自「利用有限差分法 (Finite Differential Method, FDM」· M 代码 · 共 89 行

%初始化数据; 
Nx=10;                    %划分网格  
 Ny=10;
 Lx=3;
 Ly=3;
 Hx=Lx/Nx;
 Hy=Ly/Ny;
 Nmax=160;                 %设定最大迭代次数
 u=zeros(Ny+1,Nx+1);       %初始化电位矩阵，都设为0
%根据各节点的坐标，赋边界点的电位值；
for i=1:Nx+1             %设定边界值
   u(1,i)=0; 
   u(Ny+1,i)=100;        %此处可改变100V为其它形式电位分布
end;
for i=2:Ny
   u(i,1)=0;
   u(i,Nx+1)=0;
end;
%用超松弛法迭代电势;
Error_max=0;Precision=1e-05;  %设定程序停止条件
alpha=1.55;                   %设定超松弛迭代法的α值
for k=1:Nmax                  %开始迭代
   Error_max=0;
   for i=2:Ny
    for j=2:Nx
      u1=u(i,j);
      u(i,j)=u(i,j)+ alpha*(u(i+1,j)+u(i,j+1)+u...
(i-1,j)+u(i,j-1)-4*u(i,j))/4;           %超松弛迭代法 
      u2=u(i,j);
      if abs(u2-u1)>Error_max
         Error_max=abs(u2-u1);          %比较迭代前后误差
      end
    end
   end
   if Error_max<Precision               %判断是否退出迭代
       break
   end
end
if Error_max<Precision 
  disp('迭代成功,迭代次数:') ;          %显示迭代是否成功
  disp(k) ;  disp('最大误差:');
  disp(Error_max);
else
  disp('迭代失败!') ; 
end;
%计算电场强度;
for i=2:Ny
   for j=2:Nx
     Ex(i,j)=0.5*(u(i,j-1)-u(i,j+1))/Hx;  %计算一般电场强度
     Ey(i,j)=0.5*(u(i-1,j)-u(i+1,j))/Hy;
   end;
 end;
 Ex(1,1)=-(u(1,2)-u(1,1))/Hx;             %计算边界电场强度
 Ey(1,1)=-(u(2,1)-u(1,1))/Hy;
 Ex(1,Nx+1)=-(u(1,Nx+1)-u(1,Nx))/Hx;
 Ey(1,Nx+1)=-(u(2,Nx+1)-u(1,Nx+1))/Hy;
 Ex(1+Ny,1)=-(u(1+Ny,2)-u(1+Ny,1))/Hx;
 Ey(1+Ny,1)=-(u(1+Ny,1)-u(Ny,1))/Hy;
 Ex(1+Ny,1+Nx)=-(u(1+Ny,1+Nx)-u(1+Ny,Nx))/Hx;
 Ey(1+Ny,1+Nx)=-(u(1+Ny,1+Nx)-u(Ny,1+Nx))/Hy;
 for i=1:(Ny+1)
   for j=1:(Nx+1)
      E(i,j)=sqrt(Ex(i,j)^2+Ey(i,j)^2); %计算合成场强
   end;
 end;
%绘制三维电位、电场强度分布图，二维等位线、电力线;
x0=0;y0=0;                    %设定图像起始点
for i=1:(Nx+1)                %将X坐标轴划分为Nx等份
    xu(i)=x0+(i-1)*Hx;
end;
for i=1:(Ny+1)                %将Y坐标轴划分为Ny等份
    yu(i)=y0+(i-1)*Hy;
end;
[X,Y]=meshgrid(xu,yu);        %产生以xu和xy为基准的栅格点
subplot(4,2,1);surf(X,Y,u);title('电位分布图');       
%将窗口划分为八个区域，并在第一个区域绘制电位分布图
subplot(4,2,5);surf(X,Y,abs(Ex)) ; 
%绘制x方向的电场强度布图
subplot(4,2,6);surf(X,Y,abs(Ey));   
%绘制y方向的电场强度布图
subplot(4,2,4);surf(X,Y,E); title('场强分布图');   
%绘制合场强强度布图
subplot(4,2,7);contour3(X,Y,u,10);   
%绘制三维空间中的等电位线
subplot(4,2,2);contour(X,Y,u,10); title('等电位线');    
%绘制平面等电位线
[Dx,Dy]=gradient(u,Hx,Hy);     
subplot(4,2,3);quiver(X,Y,-Dx,-Dy); title('电力线');  
%产生梯度向量并绘制平面电力线 将每一点的梯度向量以小箭头表示