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

10 views (last 30 days)
Tom Willy
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

Answers (1)

ChrisR
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 Comments
Walter Roberson
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.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by