%% Example 5.8
% Iterative Solution of an Elliptic Equation

%% Setup
close all;

N = 20;  % points in the each direction

x0 = -1;
xN = 1;
y0 = -1;
yN = 1;

dx = (xN-x0)/N;
dy = (yN-y0)/N;

x = linspace(x0,xN,N+1);
y = linspace(y0,yN,N+1);

[X Y] = meshgrid(x,y);

Q = 2*(2-X.^2 - Y.^2);
Qint = Q(2:N,2:N);  % obtain interior
q = Qint(:);        % transform into a vector

%% Direct Solve
% construct system.  
% See Applied Numerical Linear Algebra by James Demmel
% footnote on page 275
TN = gallery('tridiag',N-1,-1,2,-1);
TNxN = kron(speye(N-1),TN) + kron(TN,eye(N-1));
tic
pd = (1/(dx^2)*TNxN)\q;
toc
Pd = zeros(N+1,N+1);
Pd(2:N,2:N) = reshape(pd,N-1,N-1);

%% True solution
Pd = (X.^2-1).*(Y.^2-1);

tol = 1e-4;

%% Point Jacobi
Pj = zeros(N+1,N+1);   % initial guess
Pnew = zeros(N+1,N+1); % extra storage
err = max(max(abs(Pd(2:N,2:N)-Pj(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Pnew(i,j) = (Pj(i+1,j)+Pj(i-1,j)+Pj(i,j+1)+Pj(i,j-1))/4 + dx^2/4*Q(i,j); 
        end
    end
    Pj = Pnew;
    err = max(max(abs(Pd-Pj)))/max(max(abs(Pd)));
    Ej(iter) = err;
end

pj_iter = iter

%% Gauss-Seidel
Pgs = zeros(N+1,N+1);   % initial guess
err = max(max(abs(Pd(2:N,2:N)-Pgs(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Pgs(i,j) = (Pgs(i+1,j)+Pgs(i-1,j)+Pgs(i,j+1)+Pgs(i,j-1))/4 + dx^2/4*Q(i,j); 
        end
    end
    err = max(max(abs(Pd-Pgs)))/max(max(abs(Pd)));
    Egs(iter) = err;
end

pgs_iter = iter

%% SOR
Psor = zeros(N+1,N+1);   % initial guess
Pold = zeros(N+1,N+1);   % old solution
omega = 1.8; % SOR parameter
err = max(max(abs(Pd(2:N,2:N)-Psor(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Psor(i,j) = omega*((Psor(i+1,j)+Psor(i-1,j)+Psor(i,j+1)+Psor(i,j-1))/4 + dx^2/4*Q(i,j))+(1-omega)*Psor(i,j); 
        end
    end
    err = max(max(abs(Pd-Psor)))/max(max(abs(Pd)));
    Esor(iter) = err;
end

psor_iter = iter