% solve del^2 u = -4 * pi * rho using a fast Poisson solver
% domain: for simplicity, a square {0 <= x,y <= 1}. the edges are grounded, so
% on the boundary u(x,y) = 0. we’ll have 100 cells in the mesh of our domain.
a = 0; b = 1;
ncells = 100;
xe = linspace(a,b,ncells+1); ye = linspace(a,b,ncells+1);
[x, y] = meshgrid(xe,ye);
% grid spacing
h = (b-a)/ncells;
% set up right hand side. our charge density will just be 10 point charges on
% random grid points with strength 1. note that the right hand side does not
% include the borders; the charge is always zero on the borders since they are
% have no potential (they are grounded).
rho = zeros(ncells-1);
pts = round(rand(10,2) * (ncells-1));
for i=1:10
rho(pts(i,1),pts(i,2)) = -4 * pi;
end
% now solve to get a potential field
V = fft_poisson(rho,h);
% pad V with zero boundary values
V = [zeros(ncells-1,1), V, zeros(ncells-1,1)];
V = [zeros(1,ncells+1); V; zeros(1,ncells+1)];
% contour plot of potential field
figure(1)
contour(x,y,V)
xlabel('x'); ylabel('y');
% mesh plot of potential field
figure(2)
mesh(x,y,V)
xlabel('x'); ylabel('y');