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

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[ChrisR]

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.

[Tom Willy]

Thanks ChrisR for your help. Could you please advice what i should add.

[Walter Roberson]

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]

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.