CoCalc -- Collaborative Calculation in the Cloud
SharedProyecto_final_analisis_numerico2 / PEFD.mOpen in CoCalc

% Nombres: -Rafael An�bal Calder�n Melara  CM11099
%          -Estely Dolores M�ndez Moreno   MM09204

% Fecha de creaci�n: 15 de Noviembre de 2016
% Esta funci�n aproxima la soluci�n de la ecuaci�n de Poisson
%      f_xx(x,y) + f_yy(x,y) = f(x,y); con x en [a,b], y en [c,d]
% sujeto a las condiciones iniciales
%      u(x,y) = g(x,y) si x=a o x=b & y pertenece a [c,d]
%      u(x,y) = g(x,y) si y=c o y=d & x pertenece a [a,b]
% El algoritmo utiliza el m�todo de diferencias finitas 
% y a su vez el m�todo de Gauss Seidel.
% Notar que se utiliza la norma infinita en nuestro c�digo, 
% pero se podr�a utilizar cualquier otra norma.
% INPUT:  -- f: Funciones
%        -- a,b,c,d: Valores frontera.
%        -- n>=3, m>=3: Enteros que sirven para determinar las 
%                   dimensiones de la cuadr�cula.
%        -- tol: Tolerancia  requerida.
%        -- N: M�ximo n�mero de iteraciones.

% OUTPUT: -- Tabla con las aproximaciones para u(x,y) y los puntos 
%            x & y.
%         -- l: n�mero de iteraciones requeridas para obtener 
%              la soluci�n.

function [] = PEFD(f,a,b,c,d,m,n,tol,N)

% h y k representan el tama�o de paso para la construcci�n de los
% puntos x_i e y_j respectivamente.
h = (b-a)/n;
k = (d-c)/m;

% Se crean x & y que son vectores que contienen las coordenadas  
% para los distintos puntos de la cuadr�cula que utiliza 
% el m�todo.
x = zeros(1,n-1);
for i = 1:n-1
    x(i) = a+i*h;
end
y = zeros(1,m-1);
for j = 1:m-1
    y(j) = c+j*k;
end

% w es la matriz de ceros que luego se ir� actualizando con 
% las aproximaciones w_ji de la soluci�n u(x_i,y_j).
w = zeros(n-1,m-1);

% Se calculan lambda y mu para hacer m�s eficiente el c�digo en t�rmios de
% las operaciones que se realizan. 
lambda = (h^2)/(k^2);
mu = 2*(1+lambda);

% El bucle while se utiliza para realizar iteraciones del m�todo
% de Gauss-Seidel Optimizado
l = 1;
while l <= N 
    
    % Aqu� se actualiza el elemento 1 de la �ltima columna de w.
    z = (-(h^2)*f(x(1),y(m-1))+fung(a,y(m-1))+lambda*fung(x(1),d)+lambda*w(1,m-2)+w(2,m-1))/mu;
    NORM = norm(z-w(1,m-1),inf);
    w(1,m-1)=z;
    
    % El siguiente for actualiza los valores de w, desde el segundo
    % al pen�ltimo elemento de w, de la �ltima columna.
    for i = 2:n-2
        z = (-(h^2)*f(x(i),y(m-1))+lambda*fung(x(i),d)+w(i-1,m-1)+w(i+1,m-1)+lambda*w(i,m-2))/mu;
        if norm(w(i,m-1)-z,inf) > NORM
            NORM = norm(w(i,m-1)-z,inf);
        end
        w(i,m-1)=z;
    end
    
    % Aqu� se actualiza el �ltimo valor de la �ltima columna de w.
    z = (-(h^2)*f(x(n-1),y(m-1))+fung(b,y(m-1))+lambda*fung(x(n-1),d)+w(n-2,m-1)+lambda*w(n-1,m-2))/mu;
    if norm(w(n-1,m-1)-z,inf) > NORM
        NORM = norm(w(n-1,m-1)-z,inf);
    end
    w(n-1,m-1) = z;
    
    % El siguiente for actualiza los valores de w, desde el segundo 
    % al pen�ltimo, de la todas las fila.
    for j = 2:m-2
        
        % Aqu� actualizan los valores de w, desde el segundo al pen�ltimo, 
        % de la primera fila.
        z = (-(h^2)*f(x(1),y(j))+fung(a,y(j))+lambda*w(1,j+1)+lambda*w(1,j-1)+w(2,j))/mu;
        if norm(w(1,j)-z,inf) > NORM
            NORM = norm(w(1,j)-z,inf);
        end
        w(1,j) = z;
        
        % Este for actualiza los valores de w, desde el segundo al pen�ltimo, 
        % de las filas 2 a la n-2.
        for i = 2:n-2
            z = (-(h^2)*f(x(i),y(j))+w(i-1,j)+lambda*w(i,j+1)+w(i+1,j)+lambda*w(i,j-1))/mu;
            if norm(w(i,j)-z,inf) > NORM
                NORM = norm(w(i,j)-z,inf);
            end
            w(i,j) =z;
        end
        
        % Aqu� se actualizan los valores de w, desde el segundo al pen�ltimo,
        % de la �ltima fila
        z = (-(h^2)*f(x(n-1),y(j))+fung(b,y(j))+w(n-2,j)+lambda*w(n-1,j+1)+lambda*w(n-1,j-1))/mu;
        if norm(w(n-1,j)-z,inf) > NORM
            NORM = norm(w(n-1,j)-z,inf);
        end
        w(n-1,j) =z;
    end
    
    % Aqu� se actualiza el primer valor de la primera fila y primera columna de w. fila 1 y columna 1 de w.
    z = (-(h^2)*f(x(1),y(1))+fung(a,y(1))+lambda*fung(x(1),c)+lambda*w(1,2)+w(2,1))/mu;
    if norm(w(1,1)-z,inf) > NORM
        NORM = norm(w(1,1)-z,inf);
    end
    w(1,1) =z;
    
    % El siguiente for actualiza los valores de w, desde el segundo al pen�ltimo
    % de la primera columna.
    for i = 2:n-2
        z = (-(h^2)*f(x(i),y(1))+lambda*fung(x(i),c)+w(i-1,1)+lambda*w(i,2)+w(i+1,1))/mu;
        if norm(w(i,1)-z,inf) > NORM
            NORM = norm(w(i,1)-z,inf);
        end    
        w(i,1) =z;
    end
    
    % Aqu� se actualiza el �ltimo valor de la ultima columna de la matriz w.
    z = (-(h^2)*f(x(n-1),y(1))+fung(b,y(1))+lambda*fung(x(n-1),c)+w(n-2,1)+lambda*w(n-1,2))/mu;
    if norm(w(n-1,1)-z,inf) > NORM
        NORM = norm(w(n-1,1)-z,inf);
    end
    w(n-1,1) =z;
    
    % Este if verifica que se ha alcanzado la tolerancia deseada.
    if NORM <= tol
        break
    end
    
    % Aqu� se actualiza el contador del bucle while.
    l = l+1;
end

% Este �ltimo if env�a un mensaje de error si se ha revasado el n�mero
% de iteraciones deseado.
if l > N
error('El numero maximo de iteraciones fue excedido') 
end
I=zeros(m-1,1);
J=zeros(m-1,1);
%imprime los t�tulos de cada columna de la tabla resultante
fprintf('      i              j              xi              xj              w(i,j) \n')
for j=1:m-1 %se llena los valores para los �ndices j en la tabla resultante
    J(j)=j;
end
for i=1:n-1 %se llenan los vectores con las coordenadas de x e y
    x1=zeros(m-1,1);
    for j=1:m-1
        I(j)=i; %se llena los valores para los indices i en la tabla resultante
        x1(j)=x(i);
    end
    A= [ I  J  x1 transpose(y) transpose(w(i,:)) ]; %guarda la matriz con 
    %las aproximaciones para cada i
    %imprime los elemento de la tabla
    fprintf([repmat('%10.5f\t', 1, size(A, 2)) '\n'], A');
end
% muestra el n�mero de iteraciones
fprintf('\n Y el numero de iteraciones necesarias fue: %d ', l) 
end % termina el codigo

% Ejecuci�n del ejemplo 2 , de la secci�n 12.1 del libro de Burden
% f = @(x,y)x*exp(y);
% a=0; b=2; c=0; d=1; m=5; n=6; tol=1e-10; N=100; 
% PEFD(f,a,b,c,d,m,n,tol,N)