function [U,x,y,iter]= EllipticEqDirichletCondOptimum(a,b,c,d,h,k,tol)
% [U,x,y,iter]= EllipticEqDirichletCondOptimum(a,b,c,d,h,k,tol)
% Elliptic Equation (pde) with Dirichlet Boundary Conditions:
% Laplacian(u)(x,y)= f(x,y) in a rectangular domain D=[a, b] x [c, d]
% u(x,y) fixed in the boundary of D. 
% Finite differences, h: step in x, k: step in y (uniform mesh).
% Iterative method: SOR method (over-relax) with optimum w and tolerance tol.
% Example: [U,x,y,iter]= EllipticEqDirichletCondOptimum(0,1,0,1,0.01,0.01,1e-5);
% iter, contourf(U',10), colorbar, xlabel('x'), ylabel('y')
% figure, [X,Y]=meshgrid(x,y); surf(X,Y,U), colorbar
% xlabel('x'), ylabel('y'), zlabel('u'), shading interp

% 0.-Inizialització de valors:
maxiter= 10000;  % número màxim d'iterats

nx= (b-a)/h;  x= a:h:b;  nx1= nx+1; % punts en x
ny= (d-c)/k;  y= c:k:d;  ny1= ny+1; % punts en y
if nx1 ~= fix(nx1) | ny1 ~= fix(ny1)
    error('El número de punts en x o en y no és enter (help EllipticEqDirichletCondOptimum)')
end

alfa= k/h;  alfa2= alfa^2; beta= 2*(1+alfa2); k2= k^2;  % coeficients
% si prenem mateix pas en x i y, alfa=1, beta=4.

[g1x,g2x,g3y,g4y]= g(x,y);  % valors a la frontera
u0= (sum(g1x+g2x)/nx1+sum(g3y+g4y)/ny1)/4; % promig valors a la frontera

U= u0*ones(nx1,ny1);  % U_ij^0 aproximació inicial (matriu)
U(:,1)  = g1x';  % condicions de frontera tipus Dirichlet
U(:,ny1)= g2x';
U(1,:)  = g3y;
U(nx1,:)= g4y;

fxy= f(x,y); % source term

% Optimum over-relax factor:
rho= (cos(pi/nx)+cos(pi/ny))/2;
omega= 2/(1+sqrt(1-rho^2))
omega1= 1-omega;
omegabeta=omega/beta;

% 1.- Bucle principal (Iteració pel mètode SOR)
% Resol el sistema lineal equació per equació.

incr= tol+1;   iter= 0;
while incr>tol & iter < maxiter
    incr= 0;
    for i= 2:nx
        for j= 2:ny  % Uij^(k+1)= (1-omega)*Uij^(k) + omega*Uij^(k+1)[Gauss-Seidel]
            Uij= omega1*U(i,j)+omegabeta*(alfa2*(U(i-1,j)+U(i+1,j))+U(i,j-1)+U(i,j+1)-k2*fxy(i,j));
            incr= max(incr,abs(U(i,j)-Uij));
            U(i,j)= Uij;
        end
    end
    iter= iter+1;
end

function fxy=f(x,y)
% Source term. If fxy= 0, Laplace equation. 
[X,Y]= meshgrid(x,y); 
fxy=0.*X; % notació punt.
% fxy=zeros(nx1,ny1); % matrix.

function [g1x,g2x,g3y,g4y]=g(x,y)
% Dirichlet Bondary conditions.
% u(x,c)=g1(x), u(x,d)=g2(x)
% u(a,y)=g3(y), u(b,y)=g4(y)

g1x= 0.*x;     % notació punt.
g2x= 0.*x;
g3y= y.*(1-y);
g4y= y.*(1-y);

% (C) Jordi Ripoll 2007.