M=6;
N=6; %網格節點數6*6=36個
U1=ones(N,M); %行列二維數組
m=5,n=5; %橫縱向網格數
U1(1,:)=ones(1,M)*50; %條件邊界值
U1(N,:)=ones(1,M)*100;
for i= 1:N
U1(i,1)=0;
U1(i,M)=100;
end
t1=(cos(pi/m)+cos(pi/n))/2;
w=2/(1+sqrt(1-t1*t1));
U2=U1; P=1;T=0; %初始化
k=0
while(P>1e-5) %由v1叠代，算出v2，叠代精度1e-5
k=k+1; %計算叠代次數
P=0;
for i=2:N-1; %行循環
for j=2:M-1; %列循環
U2(i,j)=U1(i,j)+(U1(i,j+1)+U1(i+1,j)+U2(i-1,j)+U2(i,j-1)-4*U1(i,j))*w/4; %差分方程
T=abs(U2(i,j)-U1(i,j));
if (T>P) P=T; end
end
end
U1=U2;
end
subplot(1,2,1),mesh(U2); %三維圖
axis([0,6,0,6,0,100]);
subplot(1,2,2),contour(U2,15); %等電位線
hold on;
x=1:1:M; y=1:1:N
[xx,yy]=meshgrid(x,y); %柵格
[Gx,Gy]=gradient(U2,0.6,0.6); %梯度
quiver(xx,yy,Gx,Gy,-1.0,'r'); %根據梯度畫箭頭
axis([-1.5,M+2.5,-2,13]); %坐標邊框設置
plot([1,1,M,M,1],[1,N,N,1,1],'K'); %畫導體邊框
text(M/2-0.5,N+0.4,'100V','fontsize',6);%上標註
text(M/2,0.3,'50V','fontsize',6);%下標註
text(-0.3,N/2,'0V','fontsize',5);%左標註
text(M+0.1,N/2,'100V','fontsize',5);%右標註
hold off