# Getting error in code - program to solve Poisson’s equation in 1D

9 views (last 30 days)
Tom Willy on 28 Oct 2021
Edited: Walter Roberson on 5 Nov 2021
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 CommentShowHide None
Rena Berman on 5 Nov 2021
(Answers Dev) Restored edit

Sign in to comment.

### Answers (1)

ChrisR on 29 Oct 2021
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 CommentsShowHide 2 older comments
Walter Roberson on 29 Oct 2021
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.

Sign in to comment.

### Categories

Find more on Matrix Indexing in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!