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



#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define PARTICION 11  // Tamaño de la particion
#define VISUALIZA 1   // (0) No visualiza la salida, otro valor la visualiza


// Lado derecho de la ecuacion diferencial parcial
double LadoDerecho(double x)
{
   double pi = 3.1415926535897932384626433832;
   return -pi * pi * cos(pi * x);
}

// Resuelve Ax=b usando el metodo Jacobi
void Jacobi(double **A, double *x, double *b, int n, int iter)
{
   int i, j, m;
   double sum;
   double *xt;
   xt = (double*)malloc(n*sizeof(double));

   for (m = 0; m < iter; m ++)
   {
      for (i = 0; i < n; i++)
      {
         sum = 0.0;
         for (j = 0; j < n; j ++)
         {
            if ( i == j) continue;
            sum += A[i][j] * x[j];
         }
         xt[i] = (1.0 / A[i][i]) * (b[i] - sum);
      }
      for (i = 0; i < n; i++) x[i] = xt[i];
   }
   if (VISUALIZA)
   {
      printf("\nMatriz\n");
      for (i = 0; i < n; i++)
      {
         for (j = 0; j < n; j++) printf("%f ", A[i][j]);
         printf("\n");
      }
      printf("\nb\n");
      for (i = 0; i < n; i++) printf("%f ", b[i]);
      printf("\nx\n");
      for (i = 0; i < n; i++) printf("%f ", x[i]);
      printf("\n");
   }
   free(xt);
}


// Resuelve Ax=b usando el metodo Gauss-Seidel
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 ++)
   {
      for (i = 0; i < n; i++)
      {
         sum = 0.0;
         for (j = 0; j < n; j ++)
         {
            if ( i == j) continue;
            sum += A[i][j] * x[j];
         }
         x[i] = (1.0 / A[i][i]) * (b[i] - sum);
      }
   }
   if (VISUALIZA)
   {
      printf("\nMatriz\n");
      for (i = 0; i < n; i++)
      {
         for (j = 0; j < n; j++) printf("%f ", A[i][j]);
         printf("\n");
      }
      printf("\nb\n");
      for (i = 0; i < n; i++) printf("%f ", b[i]);
      printf("\nx\n");
      for (i = 0; i < n; i++) printf("%f ", x[i]);
      printf("\n");
   }
}



int main(int argc, const char* argv[])
{

   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 = PARTICION;              // Particion
   int N = n-2;                    // Incognitas
   double h = (xf - xi) / (n - 1); // Incremento en la malla
   int i;

   // Declaracion de la matriz A y los vectores b y x
   double **A;  // Matriz A
   double *b;     // Vector b
   double *x;     // Vector x
   // Solicitud de la memoria dinamica
   b = (double*)malloc(N*sizeof(double));
   x = (double*)malloc(N*sizeof(double));
   A = (double**)malloc(N*sizeof(double*));
   for (i= 0; i < N; i++) A[i] = (double*)malloc(N*sizeof(double));

   // 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][0] = P;
   A[0][1] = 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][i - 1] = R;
      A[i][i] = P;
      A[i][i + 1] = Q;
      b[i] = LadoDerecho(xi + h * (i + 1));
   }
   // Renglon final de la matriz A y vector b
   A[N - 1][N - 2] = R;
   A[N - 1][N - 1] = P;
   b[N - 1] = LadoDerecho(xi + h * N - 2) - vf * Q;

   // Resuleve el sistema lineal Ax=b por Gauss-Seidel
   for(i = 0; i < N; i++) x[i] = 0.0;
   Gauss_Seidel(A, x, b, N, N*7);
   printf("%f %f\n", xi, vi);
   for(i = 1; i < N+1; i++) printf("%f %f\n", xi + h*i, x[i-1]);
   printf("%f %f\n\n\n\n", xf, vf);


   // Resuleve el sistema lineal Ax=b por Jacobi
   for(i = 0; i < N; i++) x[i] = 0.0;
   Jacobi(A, x, b, N, N*14);
   printf("%f %f\n", xi, vi);
   for(i = 1; i < N+1; i++) printf("%f %f\n", xi + h*i, x[i-1]);
   printf("%f %f\n", xf, vf);

   // Liberacion de la memoria dinamica
   free(b);
   free(x);
   for (i= 0; i < N; i++) free(A[i]);
   free(A);

   return 0;
}