%% BVP / Elliptic PDE examples (Steady-State Heat Eqn.)
% Solving Laplace equation Txx+Tyy=0 subject to boundary conditions (9.1)

%Distance in Centimeters
w=10;       %Width
h=15;       %Height
T0=100;     %Temperature Peak
dx=2.5;    
dy=dx;
x=[0:dx:w]; 
y=[0:dy:h];
[X,Y]=meshgrid(x,y);
display('Analytical');
T=T0.*sinh(pi*Y./w).*sin(pi*X./w)/sinh(pi*h/w)

%Plot Solution
figure(1);
imagesc(x,y,T); 
axis xy;
xlabel('Width (cm)');
ylabel('Height (cm');
title('Analytical Solution: Temperature (C)');
colorbar;

%% Apply 2D PDE solution, iteratively, via Jacobi Iteration
n=length(x);
m=length(y);
Ti=zeros(size(X));
Ti(end,:)=T0*sin(pi*x/(w));
T=zeros(size(X));
toler=1e-2;
iter=0;
while abs(norm(Ti)-norm(T))>toler
    T=Ti;
    iter=iter+1;
    for i=2:m-1
        for j=2:n-1
            Ti(i,j)=T(i,j)+(1/4)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j));
        end;
    end;
end;
display('Jacobi');
Ti
iter

%Plot Results
figure(3);
imagesc(x,y,Ti); 
axis xy;
xlabel('Width (cm)');
ylabel('Height (cm');
title('Iterative Numerical Solution: Temperature (C)');
colorbar;

%% Apply 2D PDE solution, iteratively, via Gauss-Seidel
n=length(x);
m=length(y);
Ti=zeros(size(X));
Ti(end,:)=T0*sin(pi*x/(w));
T=zeros(size(X));
toler=1e-2;
iter=0;
while abs(norm(Ti)-norm(T))>toler
    T=Ti;
    iter=iter+1;
    for i=2:m-1
        for j=2:n-1
            Ti(i,j)=T(i,j)+(1/4)*(Ti(i-1,j)+T(i+1,j)+Ti(i,j-1)+T(i,j+1)-4*T(i,j));
        end;
    end;
end;
display('Gauss-Seidel');
Ti
iter

%% Apply 2D PDE solution, iteratively, via SOR with omega=1.25
n=length(x);
m=length(y);
Ti=zeros(size(X));
Ti(end,:)=T0*sin(pi*x/(w));
T=zeros(size(X));
toler=1e-2;
iter=0;
omega=1.5;
while abs(norm(Ti)-norm(T))>toler
    T=Ti;
    iter=iter+1;
    for i=2:m-1
        for j=2:n-1
            Ti(i,j)=T(i,j)+omega*(1/4)*(Ti(i-1,j)+T(i+1,j)+Ti(i,j-1)+T(i,j+1)-4*T(i,j));
        end;
    end;
end;
display('SOR');
Ti
iter
