% Ejemplo del Metodo de Diferencias Finitas para resolver la ecuacion diferencial parcial:
%    Uxx=-Pi*Pi*cos(Pix)
%        xi<=U<=xf
%        U(xi)=vi y U(xf)=vf
% con solucion analitica: 
%     cos(pix)
%
% Para conocer mas del metodo ver:
%    http://132.248.182.159/acl/MDF
%
% Autor: Antonio Carrillo Ledesma
%    http://mmc.geofisica.unam.mx/acl


% Ejecutar la funcion
% [A,b,x]=mdf1DDDBand(11);
function [A,b,x] = mdf1DDDBand(n)
    xi = 0.0;         % Inicio de dominio
    xf = 1.0;         % Fin de dominio
    vi = 1.0;         % Valor en la forntera xi
    vf = -1.0;        % Valor en la frontera xf
    N=n-2;            % Nodos interiores
    h=(xf-xi)/(n-1);  % Incremento en la malla
    A=zeros(N,3);     % Matriz A
    b=zeros(N,1);     % Vector b
    
    %Stencil
    R=1/(h^2);
    P=-2/(h^2);
    Q=1/(h^2);
    
    
    % Primer renglon de la matriz A y vector b
    A(1,2)=P;
    A(1,3)=Q;
    b(1)=LadoDerecho(xi)-vi*R;
    % Renglones intermedios de la matriz A y vector b
    for i=2:N-1 
      A(i,1)=R;
      A(i,2)=P;
      A(i,3)=Q;
      b(i)=LadoDerecho(xi+h*(i));
    end
    % Relglon final de la matriz A y vector b
    A(N,1)=R;
    A(N,2)=P;
    b(N)=LadoDerecho(xi+h*N)-vf*Q;
    
    % Resuleve el sistema lineal Ax=b
    x=Gauss(A,b,N,200);
    % Prepara la graficacion
    xx=zeros(n,1);
    zz=zeros(n,1);
    for i=1:n
      xx(i)=xi+h*(i-1);
      zz(i)=SolucionAnalitica(xx(i));
    end
    
    yy=zeros(n,1);   
    yy(1)=vi;         % Condicion inicial
    for i=1:N
      yy(i+1)=x(i);
    end
    yy(n)=vf;        % Condicion inicial
    
    % Grafica la solucion de la Ecuacion Diferencial Parcial en 1D
    plot(xx,[yy,zz]);
    %plot(xx,zz);
    
    
    % Resuleve el sistema lineal Ax=b
    x=Jacobi(A,b,N,N*14);
    % Prepara la graficacion
    xx=zeros(n,1);
    zz=zeros(n,1);
    for i=1:n
      xx(i)=xi+h*(i-1);
      zz(i)=SolucionAnalitica(xx(i));
    end
    
    yy=zeros(n,1);   
    yy(1)=vi;         % Condicion inicial
    for i=1:N
      yy(i+1)=x(i);
    end
    yy(n)=vf;        % Condicion inicial
    
    % Grafica la solucion de la Ecuacion Diferencial Parcial en 1D
    plot(xx,[yy,zz]);
    %plot(xx,zz);    
endfunction

% Lado derecho de la ecuacion
function y=LadoDerecho(x)
   y=-pi*pi*cos(pi*x);
endfunction

% Solucion analitica a la ecuacion
function y=SolucionAnalitica(x)
  y=cos(pi*x);
endfunction



% Resuelve Ax=b usando el metodo Jacobi
function x=Jacobi(A,b,N,iter)
x=zeros(N,1);
xt=zeros(N,1);
for it = 1: iter
   sum =  A(1,3) * x(2);
   xt(1) = (1.0 / A(1,2)) * (b(1) - sum);
   for i = 2: N-1
      sum = A(i,1) * x(i-1) +  A(i,3) * x(i+1);
      xt(i) = (1.0 / A(i,2)) * (b(i) - sum);
   end
   sum = A(N,1) * x(N-1);
   xt(N) = (1.0 / A(N,2)) * (b(N) - sum);
   for i = 1: N
      x(i)=xt(i);
   end
end
endfunction


% Resuelve Ax=b usando el metodo Gauss-Seidel
function x=Gauss(A,b,N,iter)
x=zeros(N,1);
for it = 1: iter
   sum =  A(1,3) * x(2);
   x(1) = (1.0 / A(1,2)) * (b(1) - sum);
   for i = 2: N-1
      sum = A(i,1) * x(i-1) +  A(i,3) * x(i+1);
      x(i) = (1.0 / A(i,2)) * (b(i) - sum);
   end
   sum = A(N,1) * x(N-1);
   x(N) = (1.0 / A(N,2)) * (b(N) - sum);
end
endfunction




% Ejecutar la funcion
% [A,b,x]=mdf1DDDBand(11);