#include "chapter13.cpp"

typedef double runtype;
typedef std::vector<runtype> Array;
typedef std::vector<Array> Matrix;

#define pi 3.1415927

//parabolic1 functions
inline runtype f(runtype x) { return (sin(x)); }


// hyperbolic functions
inline runtype f2(runtype x)
{
  return(sin(pi*x));
}

inline runtype true_solution(runtype x, runtype t)
{
  return(sin(pi*x) * cos (pi*t));
}


inline runtype f3(runtype x, runtype y) { return - 25; }
inline runtype g3(runtype x, runtype y) { return 0; }
inline runtype true_solution3(runtype x, runtype y) { return 0.5*(cosh(5*x) + cosh(5*y))/cosh(5); }
inline runtype bndy3(runtype x, runtype y)    { return true_solution3(x, y); }
inline runtype ustart3(runtype x, runtype y)  { return  1; }




ofstream outfile1("parabolic1.txt");
ofstream outfile2("parabolic2.txt");
ofstream outfile3("hyperbolic.txt");
ofstream outfile4("seidel.txt");
ofstream outfile5("elliptic.txt");


void main()
{
// parabolic1
	const int n = 10, m = 20;
	runtype h = 0.1, k = 0.005;
	const runtype u0 = 0,v0 = 0,un = 0,vn = 0;
	Array u(n+1,0); 
	Array v(n+1,0);
	int i, j;
	runtype t;

	for (i = 1; i <= n - 1; i++)
	{
		u[i]= sin(pi*i*h);
		outfile1 << "u[" << i << "] = " << u[i] << endl;
	}
	for (j = 1; j <= m; j++)
	{
		for (i = 1; i <= n - 1; i++)
		{
			v[i] = (u[i-1] + u[i+1])/2;
			outfile1 << "v[" << i << "] = " << v[i] << endl;
		}
		t = j * k;
		for (i = 1; i <= n - 1; i++)
		{
			u[i] = exp(-(pi*pi)*t) * sin(pi*i*h) - v[i];
			outfile1 << "u[" << i << "] = " << u[i] << endl;
		}
		for (i = 1; i <= n - 1; i++)
			u[i] = v[i];
	}

// parabolic2

  Array c(n,0);
  Array d(n,0);
  Array x(n,0);
  Array y(n,0);
  runtype r, s;

  s = h * h / k;
  r = 2 + s;

  for (i = 1; i <= n - 1; i++)
  {
     d[i] = r;
     c[i] = -1.0;
     x[i] = sin(pi*i*h);
	 outfile2 << "x[" << i << "] = " << x[i] << endl;
  }

  for (j = 1; j <= m; j++)
  {
    for (i = 1; i <= n - 1; i++)
     {
       d[i] = r;
       y[i] = s * x[i];
     }

     tri(n-1, c, d, c, v, v);
     for (i = 1; i <= n - 1; i++)
		 outfile2 << "y[" << i << "] = " << y[i] << endl;
     t = j * k;
     for (i = 1; i <= n - 1; i++)
     {
       x[i] = exp(-(pi*pi)*t) * sin(pi*i*h) - y[i];
		 outfile2 << "x[" << i << "] = " << x[i] << endl;

     }
     for (i = 1; i <= n - 1 ; i++)
       x[i] = y[i];
  }



// hyperbolic
  k = 0.05;
  Array w(n+1,0);
  runtype z, ro;

  u[0] = 0; v[0] = 0; w[0] = 0;  
  u[n] = 0; v[n] = 0; w[n] = 0;
  ro = (k/h) * (k/h);
  for (i = 1; i <= n - 1; i++)
    {
      z = i * h;
      w[i] = f2(z);
      v[i] = 0.5*ro*(f2(z - h) + f2(z + h)) + (1 - ro)*f2(z);
    }  

  for (j = 2; j <= m; j++)
    {
	  outfile3 << endl << endl << endl << "Numerical Solution " << endl;
      for (i = 1; i <= n - 1; i++)
	{
	  u[i] = ro*(v[i+1] + v[i-1]) + 2*(1 - ro)*v[i] - w[i];
	  outfile3 << j << ": u[" << i << "] = " << u[i] << endl;
	}
	  outfile3 << endl << "True Solution " << endl;
      for (i = 1; i <= n - 1; i++)
	{
	  w[i] = v[i];
	  v[i] = u[i];
	  t = j * k;
	  z = i * h;
	  u[i] = true_solution(z,t); // - v[i];	
	  outfile3 << j << ": u[" << i << "] = " << u[i] << endl;
	}
    }  

// seidel

  const int nx = 8,
            ny = 8,
	    itmax = 20;
  runtype ax = 0,
        bx = 1,
        ay = 0,
        by = 1;
  Matrix u3(nx+1,Array(nx+1, 0));
  runtype h3, x3, y3;

  
  h3 = (bx - ax)/nx;
  for (j = 0; j <= ny; j++)
    {
      y3 = ay + j*h3;
      u3[0][j] = bndy3(ax, y3);
      u3[nx][j] = bndy3(bx, y3);
    }
  for (i = 0; i <= nx; i++)
    {
      x3 = ax + i*h3;
      u3[i][0] = bndy3(x3, ay);
      u3[i][ny] = bndy3(x3, by);
    }

  for (j = 1; j <= ny - 1; j++)
    {
      y3 = ay + j*h3;
      for (i = 1; i <= nx - 1; i++)
	{
	  x3 = ax + i*h3;
	  u3[i][j] = ustart3(x3,y3);
	  outfile4 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;
	}
    }

  seidel(ax, ay, nx, ny, h3, itmax, u3);
  outfile4 << "itmax = " << itmax << endl;
  for (i = 1; i <= nx - 1; i++)
    for (j = 1; j <= ny - 1; j++)
	  outfile4 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;

  for (j = 0; j <= ny; j++)
    {
      y3 = ay + j*h3;
      for (i = 0; i <= nx; i++)
	{
	  x3 = ax + i*h3;
	  u3[i][j] = fabs(true_solution3(x3,y3) - u3[i][j]);
	}
    }
  outfile4 << "itmax = " << itmax << endl;
   for (i = 1; i <= nx - 1; i++)
    for (j = 1; j <= ny - 1; j++)
	  outfile4 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;

// elliptic

  ax = 0,
  bx = 1,
  ay = 0,
  by = 1;

  
  h3 = (bx - ax)/nx;
  for (j = 0; j <= ny; j++)
    {
      y3 = ay + j*h3;
      u3[0][j] = bndy3(ax, y3);
      u3[nx][j] = bndy3(bx, y3);
    }
  for (i = 0; i <= nx; i++)
    {
      x3 = ax + i*h3;
      u3[i][0] = bndy3(x3, ay);
      u3[i][ny] = bndy3(x3, by);
    }

  for (j = 1; j <= ny - 1; j++)
    {
      y3 = ay + j*h3;
      for (i = 1; i <= nx - 1; i++)
	{
	  x3 = ax + i*h3;
	  u3[i][j] = ustart3(x3,y3);
	  outfile5 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;
	}
    }

  seidel(ax, ay, nx, ny, h3, itmax, u3);
  outfile5 << "itmax = " << itmax << endl;
  for (i = 1; i <= nx - 1; i++)
    for (j = 1; j <= ny - 1; j++)
	  outfile5 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;

  for (j = 0; j <= ny; j++)
    {
      y3 = ay + j*h3;
      for (i = 0; i <= nx; i++)
	{
	  x3 = ax + i*h3;
	  u3[i][j] = fabs(true_solution3(x3,y3) - u3[i][j]);
	}
    }
  outfile5 << "itmax = " << itmax << endl;
   for (i = 1; i <= nx - 1; i++)
    for (j = 1; j <= ny - 1; j++)
	  outfile5 << "u[" << i << "][" << j << "] = " << u3[i][j] << endl;


	
}