index error problem how can fix it?

5 次查看（过去 30 天）
ljaseon
ljaseon 2021-5-15
评论： ljaseon 2021-5-15
This is the code that what i put in but there is a some error the line at abc2(a_tic)=Vexact(xind,yind,zind)
error says that Index at position 3 exceeds array boundaries (must not exceed 11).
but i have no idea to fix it
how can i do it?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FINITE DIFFERENCE SOLUTION OF POISSON'S EQUATION:
% Vxx+Vyy+Vzz=G
% USING THE METHOD OF SUCCESSIVE OVER-RELAXATION
%
% NX : NO.OF INTERVALS ALONG X-AXIS
% NY : NO.OF INTERVALS ALONG Y-AXIS
% NZ : NO.OF INTERVALS ALONG Z-AXIS
%A X B : DIMENSION OF THE SOLUTION REGION
% V(I,J,K) : POTENTIALAT GRID POINT(X,Y Z)=H*(I,J,K)
% WHERE I=0,1,......,NX,J=0,1,......,NY K=0,1,......,NZ
% H :MESH SIZE
%***************************************************************
%SPECIFY BOUNDARY VALUES TO ZEROS OFR TO AVERAGE OF FIXED VALUES
A=1; B=1; C=1;
V1=0; V2=10; V3=20; V4=-10; V5=20; V6=40;
NX=20; % 4 12 20
NY=NX;
NZ=NX;
H=A/NX;
%SET INITIAL GUESS EQUAL TO ZEROS OR TO AVERAGE OF FIXED VALUES
for I=1:NX-1
for J=1:NY-1
for K=1:NZ-1
V(I+1,J+1,K+1)=(V1+V2+V3+V4+V5+V6)/6.0;
end
end
end
%SET POTENTIALS AT FIXED NODES
for K=1:NZ-1
for J=1:NY-1
V(1,J+1,K+1) = V1;
V(NX+1,J+1,K+1)=V4;
end
end
for I=1:NX-1
for K=1:NZ-1
V(I+1,1,K+1)=V2;
V(I+1,NY+1,K+1)= V5;
end
end
for I=1:NX-1
for J=1:NY-1
V(I+1,J+1,1)=V3;
V(I+1,J+1,NZ+1)=V6;
end
end
V(1,1,1)=(V1+V2+V3)/3.0;
V(NX+1,1)=(V1+V4)/2.0;
V(1,NY+1,1)=(V2+V5)/2.0;
V(1,1,NZ+1)=(V3+V6)/2.0;
V(NX+1,NY+1,NZ+1)=V4+V5+V6/3.0;
%FIND THE OPTIMUM OVER-RELAXATION FACTOR
T=cos(pi/NX)+cos(pi/NY)+cos(pi/NZ);
W=(8-sqrt(64-16*T^2))/(T^2);
disp(['SOR Factor Omega=',num2str(W)])
W6=W/6;
%INERATION BEGINS
NCOUNT = 0;
loop= 1;
while loop == 1;
RMIN =0;
for I=1:NX-1
X=H*I;
for J=1:NY-1
Y=H*J;
for K=1:NZ-1
Z=H*K;
G=-36.0*pi*X*(Y-1.0);
R=W6*(V(I+2,J+1,K+1)+V(I,J+1,K+1)+V(I+1,J+2,K+1)+V(I+1,J,K+1)+V(I+1,J+1,K)-6.0*V(I+1,J+1,K+1)-G*H*H*H);
RMIN=RMIN+abs(R);
V(I+1,J+1,K+1)=V(I+1,J+1,K+1)+R;
end
end
end
RMIN=RMIN/(NX*NY+NZ);
if(RMIN>=0.0001)
NCOUNT=NCOUNT +1;
if(NCOUNT>100)
loop=0;
disp('SOLUTION DOES NOT CONVERGE IN 100 INTERATIONS')
end
else %Then RMIN is less than.0001 and then solution has converged
loop=0;
disp(['Solution Converges in',num2str(NCOUNT),'iterations'])
end
end
Vnum=V;
%Grab original points a through i
abc=zeros(1,9);
a_tic=1;
vec=[0:H:1];
for ii=.25:.25:.75
for jj=.25 :.25 :.75
for kk=.25:.25:.75
xind=find(vec==ii);
yind=find(vec==jj);
zind=find(vec==kk);
%disp([xind,yind,zind])
abc(a_tic)=Vnum(xind,yind,zind);
a_tic=a_tic+1;
end
end
end
%OUTPUT THE FINITIE DIFFERENCE APPROX.RESULTS
%----------------------------------------------------------------
%CALCULATE THE EXACT SOLUTION
%
%POISSON'S EQUATION WITH HOMOGENEOUS BOUNDARY CONDITIONS
%SOLVED BY SERIES EXPANSION
for I=1:NX-1
X=H*I;
for J=1:NY-1
Y=H*J;
for K=1:NZ-1
Z=H*K;
end
end
end
SUM=0;
for M=1:10 % TAKE ONLY 10 TERMS OF THE SERIES
FM=M;
for N=1:10
FN=N;
for W=1:10
FW=W;
FACTOR1=(FM*pi/A)^2 +(FN*pi/B)^2+(FW*pi/C)^2;
FACTOR2=((-1)^(M+N+W))*(-144)*A*B*C^2/(pi^2*FM*FN*FW);
FACTOR3=1+1/W*pi;
FACTOR=FACTOR2*FACTOR3/FACTOR1;
SUM=SUM+FACTOR*sin(FM*pi*X/A)*sin(FN*pi*Y/B)*sin(FW*pi*Z/C);
end
end
end
VH=SUM;
%LAPLACE'S EQUATION WITH INHOMOGENEOUS BOUNDARY CONDITIONS
%SOLVED USING THE METHOD OF SEPARATION OF VARIABLES
C1=16*V1/pi^2;
C2=16*V2/pi^2;
C3=16*V3/pi^2;
C4=16*V4/pi^2;
C5=16*V5/pi^2;
C6=16*V6/pi^2;
SUM=0;
for K=1:10 %TAKE ONLY 10 TERMS OF THE SERIES
N=2*K-1;
for M=2*K-1;
for W=2*K-1;
AN=N;,AM=M;,AW=W;
A1=sin(AN*pi*Y/B);
A2=sin(AW*pi*Z/C);
A3=sinh(AN*AW*pi*X/B*C);
A4=AN*AW*sinh(AN*AW*A*pi/B*C);
TERM1=C1*A1*A2*A3/A4;
B1=sin(AM*pi*X/A);
B2=sin(AW*pi*Z/C);
B3=sinh(AM*AW*pi*Y/A*C);
B4=AM*AW*sinh(AM*AW*B*pi/A*C);
TERM2=C2*B1*B2*B3/B4;
D1=sin(AM*pi*X/A);
D2=sin(AN*pi*Y/B);
D3=sinh(AM*AN*pi*Z/A*B);
D4=AM*AN*sinh(AM*AN*C*pi/A*B);
TERM3=C3*D1*D2*D3/D4;
E1=sin(AN*pi*Y/B);
E2=sin(AW*pi*Z/C);
E3=sinh(AN*AW*pi*(A-X)/B*C);
E4=AN*AW*sinh(AN*AW*A*pi/B*C);
TERM4=C4*E1*E2*E3/E4;
F1=sin(AM*pi*X/A);
F2=sin(AW*pi*Z/C);
F3=sinh(AM*AW*pi*(B-Y)/A*C);
F4=AM*AW*sinh(AM*AW*B*pi/A*C);
TERM5=C5*F1*F2*F3/F4;
G1=sin(AM*pi*X/A);
G2=sin(AN*pi*Y/B);
G3=sinh(AM*AN*pi*(C-Z)/A*B);
G4=AM*AN*sinh(AM*AN*C*pi/A*B);
TERM6=C6*G1*G2*G3/G4 ;
TERM=TERM1+TERM2+TERM3+TERM4+TERM5+TERM6;
SUM=SUM+TERM;
end
end
end
VI=SUM;
Vexact(I+1,J+1,K+1)=VH+VI;
%Grab original points a through i
abc2=zeros(1,9);
a_tic=1;
vec=[0:H:1];
for ii=.25:.25:.75
for jj=.25:.25:.75
for kk=.25:.25:.75
xfind=find(vec==ii);
yfind=find(vec==jj);
zfind=find(vec==kk);
%disp([xind,yind,zind])
abc2(a_tic)=Vexact(xind,yind,zind) <------------------------------------------here is the error region
a_tic=a_tic+1;
end
end
end
figure(1);
imagesc(flipud(Vnum')),
colorbar
xlabel('x'), ylabel('y'),zlabel('z')
title('Example 3.4 : Poisson PDE')
format short g
  2 个评论
G A
G A 2021-5-15
May be you have to rename xfind, yfind and zfind in the loop with error as xind,yind and zind? Or use
abc2(a_tic)=Vexact(xfind,yfind,zfind)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by