Laplace Finite Differences

laplace.m

clear all; close all

% Input Parameters
shieldxs = [0, 10, 10*cos(pi/3), 0];
shieldys = [0, 0, 10*sin(pi/3), 0];
shieldv = 10;
corexs = [4, 4, 6, 6, 4];
coreys = [5*tan(pi/6)+1, 5*tan(pi/6)-1, 5*tan(pi/6)-1, ...
    5*tan(pi/6)+1, 5*tan(pi/6)+1];
corev = -10;
%gridsize = [129, 129];
maxits = 10e6;
err = 10e-6;

exponentRange = 3:15;
sampleV = [0.5*max(shieldxs), 0.75*max(shieldys)]; %sample voltage in real units
Vs = zeros(1, length(exponentRange));

for ei=1:length(exponentRange)
    e = exponentRange(ei);
    gridsize = [2^e+1, 2^e+1];
    % Field (real units)
    realsize = [max(shieldxs), max(shieldys)];
    xs = linspace(0, realsize(1), gridsize(1));
    ys = linspace(0, realsize(2), gridsize(2));
    dx = realsize(1)/(gridsize(1)-1);
    dy = realsize(2)/(gridsize(2)-1);
    units2grid = @(x, d) x./realsize(d).*(gridsize(d)-1);

    % Grid (pixels)
    Vgrid = ones(gridsize);
    gridxs = repmat(0:gridsize(1)-1, 1, gridsize(2));
    gridys = repmat(0:gridsize(2)-1, gridsize(1), 1);
    gridys=gridys(:)';

    % Set grid initial values
    Vgrid(:) = (shieldv+corev) / 2; % set initial est V
    [in,on] = inpolygon(gridxs', gridys', ...
        units2grid(shieldxs, 1), units2grid(shieldys, 2));
    shieldpts = ~in|on;
    Vgrid(shieldpts) = shieldv; % set shield V
    [in,on] = inpolygon(gridxs, gridys, ...
        units2grid(corexs, 1), units2grid(coreys, 2));
    corepts = in;
    Vgrid(corepts) = corev; % set core V

    % iterations
    i = 2:gridsize(1)-1;
    j = 2:gridsize(2)-1;
    for it=1:maxits
        old = Vgrid;

        Vgrid(i,j) = 0.5*(dx*(Vgrid(i, j-1)+Vgrid(i, j+1)) + ...
            dy*(Vgrid(i-1, j)+Vgrid(i+1, j))) / (dx+dy);
        Vgrid(shieldpts) = shieldv; % reset shield V
        Vgrid(corepts) = corev; % reset core V

        if (max(max(abs(Vgrid-old))) < err) % convergence
            break
        end
    end
    disp(['broke after ' num2str(it) ' its']);
    Vs(ei) = Vgrid(units2grid(sampleV(1), 1), units2grid(sampleV(2), 2));
    disp(['V at (' num2str(sampleV(1)) ', ' num2str(sampleV(2)) ...
        ') is ' num2str(Vs(ei))]);

    % Calculate E
    Ex = zeros(size(Vgrid));
    Ey = zeros(size(Vgrid));
    Ex(i,j) = (Vgrid(i-1,j)-Vgrid(i+1,j)) / (2*dx);
    Ey(i,j) = (Vgrid(i,j-1)-Vgrid(i,j+1)) / (2*dy);
    Ex(shieldpts) = 0; % reset shield E
    Ey(shieldpts) = 0;
    Ex(corepts) = 0; % reset core E
    Ey(corepts) = 0;
    E = sqrt(Ex.^2 + Ey.^2);

    % Plot V and E
    down = 4;
    Ex = downsample(downsample(Ex, down)', down)'; % Downsample E for plot
    Ey = downsample(downsample(Ey, down)', down)';
    figure
    hold on
    surf(xs, ys, Vgrid','LineStyle','none');
    shading interp; colorbar
    line(shieldxs,shieldys,max(max(Vgrid))*ones(size(shieldxs)), ...
        'color','w');
    line(corexs,coreys,max(max(Vgrid))*ones(size(corexs)),'color','w');
    axis tight; axis equal
    title('V field'); xlabel('x (cm)'); ylabel('y (cm)')
    figure
    hold on
    quiver(downsample(xs, down), downsample(ys, down), Ex', Ey');
    contour(xs, ys, E', 4);
    colorbar
    plot(shieldxs, shieldys, 'color','k');
    plot(corexs, coreys, 'color','k');
    axis tight; axis equal
    title('E field: contours show lines of constant |E|');
    xlabel('x (cm)'); ylabel('y (cm)')
    
    figure(3)
    plot(10./2.^exponentRange(1:ei), Vs(1:ei), 'o-')
end