// 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



// Carga y ejecuta la funcion
// exec('./mdf1DDD.sci',-1)
// [A,b,x]=mdf1DDD(11);
function [A,b,x] = mdf1DDD(n)
    xi = 0           // Inicio de dominio
    xf = 1;          // Fin de dominio
    vi = 1;          // Valor en la forntera xi
    vf = -1;         // Valor en la frontera xf
    N=n-2;           // Nodos interiores
    h=(xf-xi)/(n-1); // Incremento en la malla
    A=zeros(N,N);    // 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,1)=P;
    A(1,2)=Q;
    b(1)=LadoDerecho(xi)-vi*R;
    // Renglones intermedios de la matriz A y vector b
    for i=2:N-1 
      A(i,i-1)=R;
      A(i,i)=P;
      A(i,i+1)=Q;
      b(i)=LadoDerecho(xi+h*i);
    end
    // Renglon final de la matriz A y vector b
    A(N,N-1)=R;
    A(N,N)=P;
    b(N)=LadoDerecho(xi+h*N)-vf*Q;
    
    // Resuleve el sistema lineal Ax=b
    x=Gauss(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
    plot2d(xx,[yy,zz]);


    // Resuleve el sistema lineal Ax=b
    x=Jacobi(A,b,N,N*7);


    // 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
    plot2d(xx,[yy,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
   for i = 1: N
      sum = 0.0;
      for j = 1: N
         if i == j then
            continue;
         end
         sum = sum + A(i,j) * x(j);
      end
      xt(i) = (b(i) - sum) / A(i,i);
   end
   sw = 0;
   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
   for i = 1: N
      sum = 0.0;
      for j = 1: N
         if i == j then
            continue;
         end
         sum = sum + A(i,j) * x(j);
      end
      x(i) = (b(i) - sum) / A(i,i);
   end
end
endfunction

// Carga y ejecuta la funcion
// exec('./mdf1DDD.sci',-1)
// [A,b,x]=mdf1DDD(11);