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


using System;
using System.IO;
using System.Text;
namespace FDM1D
{
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) {
         Console.WriteLine(" ");
         Console.WriteLine("Matriz");
         for (i = 0; i < n; i++) {
            for (j = 0; j < 3; j++) System.Console.Write(A[i,j] + " ");
            System.Console.WriteLine(" ");
         }
         System.Console.WriteLine(" ");
         System.Console.WriteLine("b");
         for (i = 0; i < n; i++) System.Console.Write(b[i] + " ");
         System.Console.WriteLine(" ");
         System.Console.WriteLine("x");
         for (i = 0; i < n; i++) System.Console.Write(x[i]+ " ");
         System.Console.WriteLine(" ");
      }
   }


   // 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) {
         Console.WriteLine(" ");
         Console.WriteLine("Matriz");
         for (i = 0; i < n; i++) {
            for (j = 0; j < 3; j++) System.Console.Write(A[i,j] + " ");
            System.Console.WriteLine(" ");
         }
         System.Console.WriteLine(" ");
         System.Console.WriteLine("b");
         for (i = 0; i < n; i++) System.Console.Write(b[i] + " ");
         System.Console.WriteLine(" ");
         System.Console.WriteLine("x");
         for (i = 0; i < n; i++) System.Console.Write(x[i] + " ");
         System.Console.WriteLine(" ");
      }
   }


   // Lado derecho de la ecuacion diferencial parcial
   public static double  LadoDerecho(double x) {
      double pi = 3.1415926535897932384626433832;
      return -pi * pi * 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,N];  // 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();
      ejem.Gauss_Seidel(A, x, b, N, 1000);
      System.Console.WriteLine(xi + " " +  vi);
      for(i = 1; i < N+1; i++) System.Console.WriteLine(xi + h*i + " " + x[i-1]);
      System.Console.WriteLine(xf + " " + vf);

      ejem.Jacobi(A, x, b, N, 1000);
      System.Console.WriteLine(xi + " " +  vi);
      for(i = 1; i < N+1; i++) System.Console.WriteLine(xi + h*i + " " + x[i-1]);
      System.Console.WriteLine(xf + " " + vf);

   }

  }
}