%  FD.m
%  本程序用矩陣實現有限差分算法
%  本程序求解第一類邊界條件下無源矩形金屬槽電位分布
%  編程人    沙威(Wei Sha) 安徽大學(Anhui University) ws108@ahu.edu.cn


tic;
clear;clc;
%%%%%%%%%%%%%%%%%
%  1.常量賦值
M=10;  %  行數
N=8;  %  列數
TOTAL=M*N;  %  所有點數
UNKNOWN=(M-2)*(N-2);  %  所有未知點數
V0=100;  %  上側電位
V1=0;    %  下側電位
V2=0;    %  左側電位 
V3=0;    %  右側電位



%%%%%%%%%%%%%%%%%%
%  2.確定系數矩陣

%%  編號從左到右，從上到下
A=zeros(TOTAL,TOTAL);  %  系數矩陣初始化

%%  上半部分
i=0;  %  系數矩陣行標號
index=N+1; %  第二行第一個元素

for m=index+1:N:TOTAL-2*N+2;  %  找到未知點每行的起始標號
    for n=m:1:m+N-3;  %  找到未知點對應的真實編號
        i=i+1;  %  系數矩陣行循環
        A(i,n)=-4;  %  中心
        A(i,n+1)=1; %  右側
        A(i,n-1)=1; %  左側
        A(i,n+N)=1; %  下側
        A(i,n-N)=1; %  上側
    end;
end;

position=i;  %  記住系數矩陣循環位置

%%  下半部分

% 1.上側
for j=1:N;  %  槽上側包括角點
    i=i+1;  %  系數矩陣行循環
    A(i,j)=1;
end;

% 1.下側
for j=TOTAL-N+1:TOTAL;  %  槽下側包括角點
    i=i+1;  %  系數矩陣行循環
    A(i,j)=1;
end;

% 1.左側
for j=N+1:N:TOTAL-2*N+1;  %  左側不包括角點
    i=i+1;  %  系數矩陣行循環
    A(i,j)=1;
end;
    
  
% 1.右側
for j=2*N:N:TOTAL-N;  %  左側不包括角點
    i=i+1;  %  系數矩陣行循環
    A(i,j)=1;
end;





%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  3.確定激勵矩陣
B=zeros(TOTAL,1);  %  激勵矩陣初始化

s=position;  %  找到矩陣需要賦值的位置

% 1.上側
for j=1:N;  %  槽上側包括角點
    s=s+1;  %  激勵矩陣行循環
    B(s,1)=V0;
end;

% 1.下側
for j=TOTAL-N+1:TOTAL;  %  槽下側包括角點
    s=s+1;  %  激勵矩陣行循環
    B(s,1)=V1;
end;

% 1.左側
for j=N+1:N:TOTAL-2*N+1;  %  左側不包括角點
    s=s+1;  %  激勵矩陣行循環
    B(s,1)=V2;
end;
    
  
% 1.右側
for j=2*N:N:TOTAL-N;  %  左側不包括角點
    s=s+1;  %  激勵矩陣行循環
    B(s,1)=V3;
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  4.  求解矩陣方程
C=A\B;

r=0;  %   初值
for p=1:M;  %  行循環
    for q=1:N  %  列循環
        r=r+1;
        VOLTAGE(p,q)=C(r,1);  %  賦值便于觀察
    end;
end;
    
disp(VOLTAGE);  %  顯示
toc;