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



public class mdf1DDDBand {
   private int Visualiza;

   // Constructor
   public mdf1DDDBand() {
      Visualiza = 1;
   }

   // Resuelve Ax=b usando el metodo Jacobi
   public void Jacobi(double[][] A, double[] x, double[] b, int n, int iter) {
      int i, j, m;
      double sum;
      double[] xt = new double [n];

      for (m = 0; m < iter; m ++) {
         sum =  A[0][2] * x[1];
         xt[0] = (1.0 / A[0][1]) * (b[0] - sum);
         for (i = 1; i < n-1; i++) {
            sum = A[i][0] * x[i-1] +  A[i][2] * x[i+1];
            xt[i] = (1.0 / A[i][1]) * (b[i] - sum);
         }
         sum = A[i][0] * x[i-1];
         xt[i] = (1.0 / A[i][1]) * (b[i] - sum);
         for (i = 0; i < n; i++) x[i] = xt[i];
      }
      if (Visualiza != 0) {
         System.out.println(" ");
         System.out.println("Matriz");
         for (i = 0; i < n; i++) {
            for (j = 0; j < 3; j++) System.out.print(A[i][j] + " ");
            System.out.println(" ");
         }
         System.out.println(" ");
         System.out.println("b");
         for (i = 0; i < n; i++) System.out.print(b[i] + " ");
         System.out.println(" ");
         System.out.println("x");
         for (i = 0; i < n; i++) System.out.print(x[i] + " ");
         System.out.println(" ");
      }
   }


   // Resuelve Ax=b usando el metodo Gauss-Seidel
   public void Gauss_Seidel(double[][] A, double[] x, double[] b, int n, int iter) {
      int i, j, m;
      double sum;

      for (m = 0; m < iter; m ++) {
         sum =  A[0][2] * x[1];
         x[0] = (1.0 / A[0][1]) * (b[0] - sum);
         for (i = 1; i < n-1; i++) {
            sum = A[i][0] * x[i-1] +  A[i][2] * x[i+1];
            x[i] = (1.0 / A[i][1]) * (b[i] - sum);
         }
         sum = A[i][0] * x[i-1];
         x[i] = (1.0 / A[i][1]) * (b[i] - sum);
      }
      if (Visualiza != 0) {
         System.out.println(" ");
         System.out.println("Matriz");
         for (i = 0; i < n; i++) {
            for (j = 0; j < 3; j++) System.out.print(A[i][j] + " ");
            System.out.println(" ");
         }
         System.out.println(" ");
         System.out.println("b");
         for (i = 0; i < n; i++) System.out.print(b[i] + " ");
         System.out.println(" ");
         System.out.println("x");
         for (i = 0; i < n; i++) System.out.print(x[i] + " ");
         System.out.println(" ");
      }
   }


   // Lado derecho de la ecuacion diferencial parcial
   public static double LadoDerecho(double x) {
      double pi = 3.1415926535897932384626433832;
      return -pi * pi * java.lang.Math.cos(pi * x);
   }


   // Funcion Principal ....
   public static void main(String[] args) {

      double xi = 0.0;               // Inicio del diminio
      double xf = 1.0;                // Fin del dominio
      double vi = 1.0;               // Valor en la frontera xi
      double vf = -1.0;                // Valor en la frontera xf
      int n = 11;                     // Particion
      int N = n - 2;                  // Nodos interiores
      double h = (xf - xi) / (n - 1); // Incremento en la malla
      int i;

      double[][] A = new double[N][3];  // Matriz A
      double[] b = new double[N];       // Vector b
      double[] x = new double[N];       // Vector x

      // Stencil
      double R = 1 / (h * h);
      double P = -2 / (h * h);
      double Q = 1 / (h * h);

      // Primer renglon de la matriz A y vector b
      A[0][1] = P;
      A[0][2] = Q;
      b[0] = LadoDerecho(xi) - vi * R;
      // Renglones intermedios de la matriz A y vector b
      for (i = 1; i < N - 1; i++) {
         A[i][0] = R;
         A[i][1] = P;
         A[i][2] = Q;
         b[i] = LadoDerecho(xi + h * (i+1));
      }
      // Renglon final de la matriz A y vector b
      A[N - 1][0] = R;
      A[N - 1][1] = P;
      b[N - 1] = LadoDerecho(xi + h * N - 2) - vf * Q;

      // Resuleve el sistema lineal Ax=b
      mdf1DDDBand ejem = new mdf1DDDBand();

      // Resuleve el sistema lineal Ax=b por Gauss-Seidel
      for(i = 0; i < N; i++) x[i] = 0.0;
      ejem.Gauss_Seidel(A, x, b, N, N*7);
      System.out.println(xi + " " +  vi);
      for(i = 1; i < N+1; i++) System.out.println(xi + h*i + " " + x[i-1]);
      System.out.println(xf + " " + vf);

      // Resuleve el sistema lineal Ax=b por Jacobi
      for(i = 0; i < N; i++) x[i] = 0.0;
      ejem.Jacobi(A, x, b, N, N*14);
      System.out.println(xi + " " +  vi);
      for(i = 1; i < N+1; i++) System.out.println(xi + h*i + " " + x[i-1]);
      System.out.println(xf + " " + vf);
  }
}