// check_nuderiv.java
// y = U(x) = x^4 +2 x^3 + 3 x^2 + 4 x + 5
// y' = Ux = 4 x^3 + 6 x^2 + 6 x + 4
// y'' = Uxx = 12 x^2 + 12 x + 6
// y''' = Uxxx = 24 x + 12
// y'''' = Uxxxx = 24
//
// consider nx specific values of x in increasing order, xg[]
// cx, cxx, cxxx, cxxxx computed by nuderiv at xg[i]  i=0..nx-1
// y' at xg[i] =    sum j=0..nx-1  cx[i][j]*U(xg[j])  first order
// y'' at xg[i] =   sum j=0..nx-1  cxx[i][j]*U(xg[j])
// y''' at xg[i] =  sum j=0..nx-1  cxxx[i][j]*U(xg[j])
// y'''' at xg[i] = sum j=0..nx-1  cxxxx[i][j]*U(xg[j])
//
// given a differential equation 
// a(x)*Uxxxx(x)  + b(x)*Uxxx(x) + c(x)*Uxx(x) + d(x)*Uxx(x) + e(x)*U(x) = f(x)// where a(x), b(x), c(x), d(x), e(x), f(x) are known and computable,
// and  U(xg[0]) and U(xg[nx-1]) boundary values are known,
// then solve for U(xg[j]) j=1..nx-2
//
// using cx, cxx, cxxx, cxxxx, a,b,c,d,e,f, U(xg[0], U(xg[5])
// build matrix A, for nx=6,  such that nx-2 equations can
// be solved for the intermediate points of U
//
// | a11 a12 a13 a14 |   | U(xg[1]) |   | f(xg[1]) |
// | a21 a22 a23 a24 | * | U(xg[2]) | = | f(xg[2]) |
// | a31 a32 a33 a34 |   | U(xg[3]) |   | f(xg[3]) |
// | a41 a42 a43 a44 |   | U(xg[4]) |   | f(xg[4]) |
//
// the first equation will be based on i=1 with 
//    a11 = sum  a(xg[1])*cx[1][1] + b(xg[1])*cxx[1][1] + c ...
//    a12 = sum  a(xg[2])*cx[1][2] + b(xg[2])*cxx[1][2] + c ...
//    a13 = sum  a(xg[3])*cx[1][3] + b(xg[3])*cxx[1][3] + c ...
//    a14 = sum  a(xg[4])*cx[1][4] + b(xg[4])*cxx[1][4] + c ...
//
//   and the right hand side becomes
//   f(xg[1]) - a(xg[0])*cxxxx[1][0] - a(xg[5])*cxxxx[1][5]
//            - b(xg[0])*cxxx[1][0]  - b(xg[5])*cxxx[1][5]
//            ...
// then do rest of equations
//
// thus, solving the simultaneous equations gives the solution for
// unknown values of U(xg[1]), U(xg[2]), U(xg[3]), U(xg[4])
// other values may then be computed by interpolation or
// fitting the solution to a polynomial.
//

// uses deriv.java  simeq.java
class check_nuderiv
{
  int nx = 6; // number of points
  int neqn = nx-2; // number of equations
  double xg[] = { -2.0, -0.5, 0.0, 0.5, 1.0, 2.0};
  double cx[][] = new double[nx][nx]; // first derivative coef at i,j
  double cxx[][] = new double[nx][nx]; // second derivative
  double cxxx[][] = new double[nx][nx]; // third derivative
  double cxxxx[][] = new double[nx][nx]; // fourth derivative

  // functions for U and derivatives
  double U(double x)
  {
    return x*x*x*x + 2.0*x*x*x + 3.0*x*x + 4.0*x + 5.0;
  } // end U

  double Ux(double x)
  {
    return  4.0*x*x*x + 6.0*x*x + 6.0*x + 4.0;
  } // end Ux

  double Uxx(double x)
  {
    return 12.0*x*x + 12.0*x + 6.0;
  } // end Uxx

  double Uxxx(double x)
  {
    return  24.0*x + 12.0;
  } // end Uxxx

  double Uxxxx(double x)
  {
    return  24.0;
  } // end Uxxxx

  // arbitrary functions for coefficients of differential equation
  double a(double x)
  {
    return 1.0 + x;
  }

  double b(double x)
  {
    return 2.0 - x;
  }

  double c(double x)
  {
    return 1.0 + x*x;
  }
  double d(double x)
  {
    return 3.0 - x*x;
  }

  double e(double x)
  {
    return 4.0 + 3.0*x + 2.0*x*x + x*x*x;
  }
  //            7      6      5       4       3       2
  //   f(x)=   x  + 4 x  + 6 x  + 26 x  + 48 x  + 42 x  + 121 x + 86
  double f(double x) // from check_nuderiv_mws.out
  {
    return Math.pow(x,7) + 4.0*Math.pow(x,6) + 6.0*Math.pow(x,5) +
           26.0*Math.pow(x,4) + 48.0*Math.pow(x,3) + 42.0*Math.pow(x,2) + 
           121*x + 86;
  }

  int S(int i) // subscripts start at zero
  {
    return i-1;
  } 

  check_nuderiv()
  {
    double cderiv, uderiv, err, maxerr;
    int nord;
    double x, y, eqn;
    double A[][] = new double[neqn][neqn];
    double F[] = new double[neqn];
    double Uc[] = new double[neqn]; // A * Uc = F
    int row, col;

    System.out.println("check_nuderiv.java in deriv.c running ");
    System.out.println("compare nuderiv to Ux, Uxx, Uxxx ");
    System.out.println("xg="+xg[0]+", "+xg[1]+", "+xg[2]+", "+
                        xg[3]+", "+xg[4]+", "+xg[5]);
    System.out.println("U="+U(xg[0])+", "+U(xg[1])+", "+U(xg[2])+", "+
                        U(xg[3])+", "+U(xg[4])+", "+U(xg[5]));
    System.out.println("Ux="+Ux(xg[0])+", "+Ux(xg[1])+", "+Ux(xg[2])+", "+
                        Ux(xg[3])+", "+Ux(xg[4])+", "+Ux(xg[5]));
    System.out.println("Uxx="+Uxx(xg[0])+", "+Uxx(xg[1])+", "+Uxx(xg[2])+", "+
                        Uxx(xg[3])+", "+Uxx(xg[4])+", "+Uxx(xg[5]));
    System.out.println("Uxxx="+Uxxx(xg[0])+", "+Uxxx(xg[1])+", "+
                        Uxxx(xg[2])+", "+Uxxx(xg[3])+", "+Uxxx(xg[4])+", "+
                        Uxxx(xg[5]));
    System.out.println("Uxxxx="+Uxxxx(xg[0])+", "+Uxxxx(xg[1])+", "+
                        Uxxxx(xg[2])+", "+Uxxxx(xg[3])+", "+Uxxxx(xg[4])+", "+
                        Uxxxx(xg[5]));

    maxerr = 0.0;
    for(int i=0; i<nx; i++)
    {
      nord = 1;
      new nuderiv(nord, nx, i, xg, cx[i]);
      System.out.println("nx="+nx+", point i="+i+", norder="+nord);
      System.out.println("pt i="+i+" cx="+cx[i][0]+", "+cx[i][1]+", "+
                         cx[i][2]+", "+cx[i][3]+", "+cx[i][4]+", "+cx[i][5]);
      uderiv = Ux(xg[i]);
      cderiv = 0.0;
      for(int j=0; j<nx; j++) cderiv += cx[i][j]*U(xg[j]);
      err = Math.abs(uderiv-cderiv);
      maxerr = Math.max(maxerr, err);
      System.out.println("Ux="+uderiv+", ceriv="+cderiv+", err="+err);

      nord = 2;
      System.out.println("nx="+nx+", point i="+i+", norder="+nord);
      new nuderiv(nord, nx, i, xg, cxx[i]);
      System.out.println("pt i="+i+" cxx="+cxx[i][0]+", "+cxx[i][1]+", "+
                         cxx[i][2]+", "+cxx[i][3]+", "+cxx[i][4]+", "+
                         cxx[i][5]);
      uderiv = Uxx(xg[i]);
      cderiv = 0.0;
      for(int j=0; j<nx; j++) cderiv += cxx[i][j]*U(xg[j]);
      err = Math.abs(uderiv-cderiv);
      maxerr = Math.max(maxerr, err);
      System.out.println("Uxx="+uderiv+", ceriv="+cderiv+", err="+err);

      nord = 3;
      System.out.println("nx="+nx+", point i="+i+", norder="+nord);
      new nuderiv(nord, nx, i, xg, cxxx[i]);
      System.out.println("pt i="+i+" cxxx="+cxxx[i][0]+", "+cxxx[i][1]+", "+
                         cxxx[i][2]+", "+cxxx[i][3]+", "+cxxx[i][4]+", "+
                         cxxx[i][5]);
      uderiv = Uxxx(xg[i]);
      cderiv = 0.0;
      for(int j=0; j<nx; j++) cderiv += cxxx[i][j]*U(xg[j]);
      err = Math.abs(uderiv-cderiv);
      maxerr = Math.max(maxerr, err);
      System.out.println("Uxxx="+uderiv+", ceriv="+cderiv+", err="+err);

      nord = 4;
      System.out.println("nx="+nx+", point i="+i+", norder="+nord);
      new nuderiv(nord, nx, i, xg, cxxxx[i]);
      System.out.println("pt i="+i+" cxxxx="+cxxxx[i][0]+", "+cxxxx[i][1]+", "+
                         cxxxx[i][2]+", "+cxxxx[i][3]+", "+
                         cxxxx[i][4]+", "+cxxxx[i][5]);
      uderiv = Uxxxx(xg[i]);
      cderiv = 0.0;
      for(int j=0; j<nx; j++) cderiv += cxxxx[i][j]*U(xg[j]);
      err = Math.abs(uderiv-cderiv);
      maxerr = Math.max(maxerr, err);
      System.out.println("Uxxxx="+uderiv+", ceriv="+cderiv+", err="+err);
    }
    System.out.println("check_nuderiv.java finished checking cx, cxx.., maxerr="+maxerr);

    // now set up differential equation
    System.out.println("");
    System.out.println("solve a(x)*Uxxxx(x)  + b(x)*Uxxx(x) + c(x)*Uxx(x) + d(x)*Ux(x) + e(x)*U(x) = f(x) ");
    System.out.println("    x      a(x)    b(x)    c(x)    d(x)    e(x)     f(x) are: ");
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      System.out.println(x+", "+a(x)+", "+b(x)+", "+c(x)+", "+d(x)+", "+
                         e(x)+", "+f(x));
    }
    System.out.println("");
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      eqn =  a(x)*Uxxxx(x)+b(x)*Uxxx(x)+c(x)*Uxx(x)+d(x)*Ux(x)+e(x)*U(x);
      err = f(x) - eqn;
      System.out.println("f("+x+")="+f(x)+", eqn="+eqn+", err="+err);
    }
    System.out.println("");

    // set up A * Uc = F  method 1
    for(int i=1; i<nx-1; i++) // equation number, subscripts start at 0
    { 
      row = S(i);
      F[row] = 0.0;
      for(int j=1; j<nx-1; j++) // col
      {
        col = S(j);
        A[row][col] = 0.0;
      } // j
    } //  i
    System.out.println("internal cells zeroed ");
    for(int i=1; i<nx-1; i++) // equation number, subscripts start at 0
    { 
      row = S(i);
      for(int j=1; j<nx-1; j++) // col
      {
        col = S(j);
        // make entries for each term in differential equation
        A[row][col] = a(xg[i])*cxxxx[i][j] + b(xg[i])*cxxx[i][j] +
	  c(xg[i])*cxx[i][j] + d(xg[i])*cx[i][j];
      } // j
      A[row][row] += e(xg[i]);
      F[row] = f(xg[i]);
      F[row] -= a(xg[i])*cxxxx[i][0]*U(xg[0]) + b(xg[i])*cxxx[i][0]*U(xg[0]) +
                c(xg[i])*cxx[i][0]*U(xg[0])   + d(xg[i])*cx[i][0]*U(xg[0]);
      F[row] -= a(xg[i])*cxxxx[i][5]*U(xg[5]) + b(xg[i])*cxxx[i][5]*U(xg[5]) +
                c(xg[i])*cxx[i][5]*U(xg[5])   + d(xg[i])*cx[i][5]*U(xg[5]);

    } //  i
    // print A and F
    System.out.println(" A and F method 1 ");
    for(int i=0; i<4; i++)
    {
      System.out.println("|"+A[i][0]+" "+A[i][1]+" "+A[i][2]+" "+A[i][3]+
                         "| |"+ F[i]+"|");
    } // i
    new simeq(nx-2, A, F, Uc);
    System.out.println("U[1]="+Uc[0]+", "+Uc[1]+", "+Uc[2]+", "+Uc[3]);
    System.out.println("U    "+U(xg[1])+", "+U(xg[2])+", "+U(xg[3])+", "+
                        U(xg[4]));
    System.out.println("");
    System.out.println("");


    // set up A * Uc = F  method 2
    for(int i=1; i<nx-1; i++) // equation number, subscripts start at 0
    { 
      row = S(i);
      F[row] = 0.0;
      for(int j=1; j<nx-1; j++) // col
      {
        col = S(j);
        A[row][col] = 0.0;
      } // j
    } //  i
    System.out.println("internal cells zeroed ");
    for(int i=1; i<nx-1; i++) // inside boundary 
    {
      row = S(i); // equation index 
      // for each term 

      // a(x)*Uxxxx 
      for(int k=0; k<nx; k++)
      {
        if(k==0 || k==nx-1)
        {
          F[row] = F[row] - a(xg[i])*cxxxx[i][k]*U(xg[k]);
        }	 
        else
        {
          A[row][S(k)] = A[row][S(k)] + a(xg[i])*cxxxx[i][k];
        }
      } // k

      // b(x)*Uxxx 
      for(int k=0; k<nx; k++)
      {
        if(k==0 || k==nx-1)
        {
          F[row] = F[row] - b(xg[i])*cxxx[i][k]*U(xg[k]);
        }	 
        else
        {
          A[row][S(k)] = A[row][S(k)] + b(xg[i])*cxxx[i][k];
        }
      } // k

      // c(x)*Uxx 
      for(int k=0; k<nx; k++)
      {
        if(k==0 || k==nx-1)
        {
          F[row] = F[row] - c(xg[i])*cxx[i][k]*U(xg[k]);
        }	 
        else
        {
          A[row][S(k)] = A[row][S(k)] + c(xg[i])*cxx[i][k];
        }
      } // k

      // d(x)*Ux
      for(int k=0; k<nx; k++)
      {
        if(k==0 || k==nx-1)
        {
          F[row] = F[row] - d(xg[i])*cx[i][k]*U(xg[k]);
        }
        else
        {
          A[row][S(k)] = A[row][S(k)] + d(xg[i])*cx[i][k];
        }
      } // k

      // e(x)*U(x) 
      A[row][row] = A[row][row] + e(xg[i]); // row=S(i) 

      // RHS constant 
      F[row] = F[row] + f(xg[i]);
      // end terms 
    } // i 

    // print A and F
    System.out.println(" A and F method 2");
    for(int i=0; i<4; i++)
    {
      System.out.println("|"+A[i][0]+" "+A[i][1]+" "+A[i][2]+" "+A[i][3]+
                         "| |"+F[i]+"|");
    } // i
    new simeq(nx-2, A, F, Uc);
    System.out.println("U[1]="+Uc[0]+", "+Uc[1]+", "+Uc[2]+", "+Uc[3]);
    System.out.println("U    "+U(xg[1])+", "+U(xg[2])+", "+U(xg[3])+", "+
                       U(xg[4]));
    System.out.println("");

    System.out.println("end check_nuderiv.java ");
  }

  public static void main (String[] args)
  {
    new check_nuderiv();
  }
} // end check_nuderiv.java