Getting error in code - program to solve Poisson’s equation in 1D
显示 更早的评论
Trying to code to solve Poisson’s equation in 1D using Equation:
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
0 个投票
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 个评论
Tom Willy
2021-10-29
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2021-11-5
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
