Getting error in code - program to solve Poisson’s equation in 1D
2 次查看(过去 30 天)
显示 更早的评论
Trying to code to solve Poisson’s equation in 1D using Equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/791194/image.png)
But keep getting error code:
Index exceeds the number of array elements. Index must not exceed 550.
Error in Exercise_Capacitance (line 110)
phiN + (rho(i+1)-2*rho(i)/rho(i-1))/6;
Kindly advise. Thanks
The code is:
% Parameters
%
% Units:
%
% length Angstroms
% energy eV
% charge e
%
% Simulation Parameters
%
% 1 F/m = 1 C/Vm = 1e-10/1.602176634e-19 e/VA = 1e9/1.602176634 e/VA = 0.624150907e9 e/VA
% e0 = 8.8541878176e-12 F/m
%
a = 0.2; % The mesh spacing (A)
e0 = 5.52634936105e-3; % The permitivitty of free space (e/VA)
%%%%%%%%%%%%%%%%%%%%%%%
% Initialize variables
%%%%%%%%%%%%%%%%%%%%%%%
% Capacitor parameters
%
C_N = 2; % The number of capacitors
C_sigma = 1.0; % Charge per unit area (e/A^2)
C_a = 8.0; % Capacitor plate width, less the charge (A)
C_b = 2.0; % The width of the charge region (A)
C_w = [5.0, 5.0]; % The spacing between the plates within a capacitor (A)
C_l = 20.0; % The spacing between capacitors (A)
% Compute the number of mesh points needed from the dimensions of the capacitors
%
N_a = round(C_a/a);
N_b = round(C_b/a);
N_w = [round(C_w(1)/a), round(C_w(2)/a)];
N_l = round(C_l/a);
N = N_l;
for c = 1:C_N
N = N + 2*N_a + 2*N_b + N_w(c) + N_l;
end
% Initialise charge density
%
rho0 = 0.0; % This needs to be zero
rhoNp1 = 0.0; % This needs to be zero
rho = zeros(N,1);
% Build charge density for each capacitor
%
rho_p = C_sigma/C_b;
p0 = 0;
for c = 1:C_N
p_left = p0 + N_l + N_a;
p_right = p_left + N_b + N_w(c);
for p = 1:N_b
rho(p_left + p) = +rho_p;
rho(p_right + p) = -rho_p;
end
p0 = p0 + N_l + 2*(N_a + N_b) + N_w(c);
end
% Initialise electrostatic field
%
E0 = 0.0;
phi0 = 0.0;
phi_Int = zeros(N,1);
% Grid
%
X = a:a:N*a;
%%%%%%%%%%%%%%%%%%%%%%%
% Integral solution
%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:N
sum = 0.0;
for j = 1:N
sum = sum + rho(j)*n2(a, X(i)-X(j));
end
phi_Int(i) = phi0 - E0*X(i) - sum/e0;
end
phiN = phi_Int(N);
%%%%%%%%%%%%%%%%%%%%%%%
% FE solution
%%%%%%%%%%%%%%%%%%%%%%%
% Build matrix for left hand side
%
M = zeros(N);
lhs = zeros(N,1);
for i = 1 : N
eqn = eqn + 1;
if ( i == 1 || i == N )
M(i,i) = 1.0;
lhs(i) = phi0;
else
M(i,i) = 2.0 / a / a;
M(i,i+1) = - 1.0 / a / a;
M(i+1,i) = - 1.0 / a / a;
lhs(i) = phi0;
end
end
% Build vector for right hand side
%
B = -rho/e0;
phi_Int(1) + phi0;
phiN + (rho(i+1)-2*rho(i)/rho(i-1))/6;
% Solve simultaneous equations: x = M^-1 B
%
phi_FE = M\B;
%%%%%%%%%%%%%%%%%%%%%%%
% Plot solutions
%%%%%%%%%%%%%%%%%%%%%%%
plot(X, phi_Int, X, phi_FE);
xlabel("x (A)");
ylabel("V (V)");
legend("Integral", "Finite element");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integral used to compute potential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function t = n2(a, x)
if x < -a
t = 0.0;
elseif x < 0.0
t = (x+a)^3/(6.0*a);
elseif x < a
t = a*x + (a-x)^3/(6.0*a);
else
t = a*x;
end
end
回答(1 个)
ChrisR
2021-10-29
Tom,
The index i runs from 1 to N, but if rho has N elements, then rho(i+1) will give an error if i = N.
3 个评论
Walter Roberson
2021-10-29
You have
for i = 1 : N
At the end of the loop, since you have no break statements, you can be sure that i will be left at the the last value, N.
The line giving you trouble is after the loop. You index by i, but i has its final loop value and so will be N. It would be clearer in your program if you were to replace the references to i in the line to instead be references to N. N-1 and N+1 instead of i-1 and i+1. If using N+1 does not fit with the sizes of the array then that means that you are using the wrong calculations.
Walter Roberson
2021-10-29
Permanent fix: comment out or delete the two lines
phi_Int(1) + phi0;
phiN + (rho(i+1)-2*rho(i)/rho(i-1))/6;
Both lines do calculations and throw away the results and so the only point in retaining the lines would be if you are trying to deliberately provoke error messages.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!