아래 코드로 진행했을때 입력 상수 문제라고 나오는데 어떤 부분이 잘 못된건가요?
function V = Poisson_FDM_Solver_2D(V,BC,EPS,RHO,h)
% Permittivity of free-space (F/m).
EPS_0 = 8.854e-12;
% Extract simulation domain size.
[Ny,Nx] = size(V,EPS_0);
% Total number of unknowns to solve for.
L = Nx*Ny;
% Instantiate the dielectric coefficient matrices.
a0 = zeros(Ny,Nx);
a1 = zeros(Ny,Nx);
a2 = zeros(Ny,Nx);
a3 = zeros(Ny,Nx);
a4 = zeros(Ny,Nx);
% Short-hand index notation.
X1 = 2:Nx-1;
Y1 = 2:Ny-1;
X2 = 1:Nx-2;
Y2 = 1:Ny-2;
% Compute the dielectric coefficients.
a0(Y1,X1) = -( EPS(Y1,X1) + EPS(Y2,X1) + EPS(Y1,X2) + EPS(Y2,X2) );
a1(Y1,X1) = (0.5)*( EPS(Y1,X1) + EPS(Y2,X1) );
a2(Y1,X1) = (0.5)*( EPS(Y1,X2) + EPS(Y1,X1) );
a3(Y1,X1) = (0.5)*( EPS(Y2,X2) + EPS(Y1,X2) );
a4(Y1,X1) = (0.5)*( EPS(Y2,X1) + EPS(Y2,X2) );
% Separate coefficients into real and imaginary components. Matlab runs a
% lot quicker if we can limit ourselves to purely real-valued variable
% manipulation until the last moment.
a0r = real(a0); % Real
a1r = real(a1); %
a2r = real(a2); %
a3r = real(a3); %
a4r = real(a4); %
a0i = imag(a0); % Imaginary
a1i = imag(a1); %
a2i = imag(a2); %
a3i = imag(a3); %
a4i = imag(a4); %
% Check for charge densities.
if nargin == 3
b = zeros(Ny*Nx,1);
else
% Insert default h-value if not specified.
if nargin == 4
h = 1;
end
% Normalized charge matrix. This term comes from integrating the
% charge density over a square region with size h.
b = -RHO*h^2/EPS_0;
b = b(:); % Vectorize.
end
% Set the four corners to Dirichlet boundaries to avoid confusion in the
% algorithm. These points do not really matter anyway.
BC(1,1) = 1;
BC(1,Nx) = 1;
BC(Ny,Nx) = 1;
BC(Ny,1) = 1;
% Define flags for Neumann boundaries. Each case needs to be handled in
% its own unique way, so it helps to give each one a unique marker.
TOP_FLAG = -1;
BOTTOM_FLAG = -2;
LEFT_FLAG = -3;
RIGHT_FLAG = -4;
% Find the Neumann BCs on the edges.
Neumann_TOP = BC(1,:) ~= 1;
Neumann_BOTTOM = BC(Ny,:) ~= 1;
Neumann_LEFT = BC(:,1) ~= 1;
Neumann_RIGHT = BC(:,Nx) ~= 1;
% Set flags for Neumann boundaries.
BC(1,Neumann_TOP) = TOP_FLAG;
BC(Ny,Neumann_BOTTOM) = BOTTOM_FLAG;
BC(Neumann_LEFT,1) = LEFT_FLAG;
BC(Neumann_RIGHT,Nx) = RIGHT_FLAG;
% Initialize indices and values to fill the system matrix.
NZmax = 5*L; % Maximum possible number of nonzero elements.
I = zeros(NZmax,1); % i-indices.
J = zeros(NZmax,1); % j-indices.
Sr = zeros(NZmax,1); % Real Values.
Si = zeros(NZmax,1); % Imag Values.
idx = 1; % Nonzero entry index.
% Begin iterating over the unknowns and prepare to fill the system
% matrix. Filling the A-matrix is MUCH faster if you know the indices and
% values "a priori" for ALL non-zero elements. So rather than directly
% fill in the A-matrix, this loop stores all the nonzero (i,j) indices
% along with their values. The A-matrix is then instantiated directly
% from this information.
% Note that by convention, Matlab scans matrices columnwise when only
% a single element is indexed. So rather than waste time vectorizing the
% BC-matrix or initial-valued V-matrix, just use a single index to load
% and store values and remember this convention.
disp('Filling System Matrix');
for n = 1:L
% Check boundary condition.
switch BC(n)
% Dirichlet boundary.
case 1
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% Specify right-hand side as the voltage at this point.
b(n) = V(n);
% Top Neumann boundary.
case TOP_FLAG
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n+1) = -1.
I(idx) = n;
J(idx) = n + 1;
Sr(idx) = -1;
idx = idx + 1;
% Bottom Neumann boundary.
case BOTTOM_FLAG
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n-1) = -1.
I(idx) = n;
J(idx) = n - 1;
Sr(idx) = -1;
idx = idx + 1;
% Left Neumann boundary.
case LEFT_FLAG
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n+Ny) = -1.
I(idx) = n;
J(idx) = n + Ny;
Sr(idx) = -1;
idx = idx + 1;
% Right Neumann boundary.
case RIGHT_FLAG
% Set A(n,n) = 1
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n-Ny) = -1.
I(idx) = n;
J(idx) = n - Ny;
Sr(idx) = -1;
idx = idx + 1;
% Regular point. Apply the 5-point star.
otherwise
% Convention for 5-point star:
%
% V2 | Indexing direction
% V3 V0 V1 |
% V4 \-/
%
% REMINDER: Single-valued indexing of a Matrix will scan
% COLUMN-WISE!
% V0 (center) term: A(n,n) = a0(n).
I(idx) = n;
J(idx) = n;
Sr(idx) = a0r(n); % Real
Si(idx) = a0i(n); % imaginary
idx = idx + 1;
% V1 (right) term: A(n,n+Ny) = a1(n)
I(idx) = n;
J(idx) = n + Ny;
Sr(idx) = a1r(n); % Real
Si(idx) = a1i(n); % imaginary
idx = idx + 1;
% V2 (top) term: A(n,n+1) = a2(n)
I(idx) = n;
J(idx) = n + 1;
Sr(idx) = a2r(n); % Real
Si(idx) = a2i(n); % imaginary
idx = idx + 1;
% V3 (left) term: A(n,n-Ny) = a3(n)
I(idx) = n;
J(idx) = n - Ny;
Sr(idx) = a3r(n); % Real
Si(idx) = a3i(n); % imaginary
idx = idx + 1;
% V4 (bottom) term: A(n,n-1) = a4(n)
I(idx) = n;
J(idx) = n - 1;
Sr(idx) = a4r(n); % Real
Si(idx) = a4i(n); % imaginary
idx = idx + 1;
end
end
% Throw out the leftover zeros.
I = I(1:idx-1);
J = J(1:idx-1);
Sr = Sr(1:idx-1);
Si = Si(1:idx-1);
% Combine real and imaginary coefficients into one vector.
S = Sr + 1j*Si;
% Fill the system matrix.
A = sparse(I,J,S,L,L,NZmax);
% Clear out memory before performing inversion. The matrix inversion
% process is HUGELY memory intensive, and will greedily hog up all the RAM
% it can get. Clearing out the major variables from the workspace will at
% least give a few extra MB of memory to this process. It won't be much,
% but every little bit helps if we're inverting a very large system.
clear I J S BC a0 a1 a2 a3 a4 a0r a0i Sr Si ;
clear a1r a1i a2r a2i a3r a3i a4r a4i Neumann_TOP Neumann_BOTTOM;
clear Neumann_LEFT Neumann_RIGHT EPS RHO;
% Invert the matrix. The slash operator is fantastic at doing this
% quickly for sparse matrices.
disp('Solving System Matrix');
V = A\b;
disp('Done!');
% Put the voltages back into a matrix.
V = reshape(V,Ny,Nx);
return;
