// pde_nl21.java  non linear PDE, second order, second degree            
//                                                                   
// The PDE to be solved for u(x) is:        
// Uxx(x) + 2*Ux(x) + U(x)^2  = f(x)
//
// for testing code:
//        f(x) = x^4 + 2x^2 + 4x + 3
//        u(xmin=0)=1  u(xmax=1)=2
// analytic solution and boundary  ub(x) = x^2+1
//        thus  Ux(x) = 2*x
//              Uxx(x) = 2 
//
//  needs  simeq_newton5.java
//  needs  invert.java
//  needs  rderiv.java

class pde_nl21
{
  // number of points for descretization         
  int nx = 5;
  int dof = (nx-2); // DOF, not including nonlinear terms
  double hx;
  double xg[] = new double[nx];
  double ug[] = new double[nx]; // computed solution, incl boundary
  double uA[] = new double[nx]; // actual solution, incl boundary

  // Define the region for the solution:                               
  // xmin <= x <= xmax,  using  nx  points      

  double xmin =  0.0;
  double xmax =  1.0;

  // coefficients of derivatives
  double cx[][] = new double[nx][nx];
  double cxx[][] = new double[nx][nx];

  simeq_newton5 SN = new simeq_newton5(1.0e-4, 1.0, 12, false);
  // can adjust parameters             eps     b    itr debug

  pde_nl21() // construct and run
  {
    double x;

    System.out.println("The PDE to be solved is:");
    System.out.println(" Uxx(x) + 2*Ux(x) + U(x)^2 = f(x)");
    System.out.println(" f(x) = x^4 + 2x^2 + 4x + 3");
    System.out.println(" discrete u(xmin=0)=1  u(xmax=1)=2");
    System.out.println(" analytic boundary ub(x) = x^2+1");
    System.out.println(" ");

    // initialize hx, xg
    hx = (xmax-xmin)/(double)(nx-1);
    for(int i=0; i<nx; i++) xg[i] = xmin+hx*(double)i;

    // initialize ug with boundary values, uA all values
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      uA[i] = ub(x); // for testing
      if(i==0 || i==nx-1 ) ug[i] = ub(x); // known boundary
      else                 ug[i] = 0.0;
      // debug
      System.out.println("xg["+i+"]="+xg[i]+", uA["+i+"]="+uA[i]);
    }
    System.out.println(" ");

    // compute coefficients of derivatives
    for(int i=0; i<nx; i++)
    {
      // compute derivative coefficients in x direction
      new rderiv(1,nx,i,hx,cx[i]);    // Ux
      new rderiv(2,nx,i,hx,cxx[i]);   // Uxx
    }

    solve_pde();

    check_soln();
  }
 
  // must compute boundary.
  // also solution for checking 
  double ub(double x)
  {
    return x*x + 1.0;
  }

  double ubx(double x)
  {
    return  2.0*x;
  }

  double ubxx(double x)
  {
    return  2.0;
  }


  // RHS, this becomes defined by PDE and ub(x)  

  double f(double x)
  {
    return  x*x*x*x + 2.0*x*x + 4.0*x + 3.0;
  } // end f

  void solve_pde()
  {
    int nlin = 3; // number of nonlinear terms after linear terms
    int nxl = dof+nlin;
    double A[][] = new double[dof][nxl];
    double Y[] = new double[dof]; // RHS
    double us[] = new double[nxl]; // initial guess solution, offset from ug
    int var1[] = new int[nxl]; // definition of nonlinear system of equations
    int var2[] = new int[nxl];
    int var3[] = new int[nxl];
    int vari1[] = new int[nxl];
    int vari2[] = new int[nxl];

    double x;
    int ia;

    for(int i=0; i<dof; i++)
    {
      var1[i] = i; // linear terms, required
      var2[i] = -1;
      var3[i] = -1;
      vari1[i] = -1;
      vari2[i] = -1;
    }
    for(int i=dof; i<nxl; i++)
    {
      // set var1, var2 for nonlinear terms
      var1[i] = i-dof; // x[0]^2 x[1]^2 x[2]^2
      var2[i] = i-dof;
      var3[i] = -1;
      vari1[i] = -1;
      vari2[i] = -1;
    }
    

    // clear A and Y
    for(int i=0; i<dof; i++)
    {
      for(int j=0; j<nxl; j++)
      {
	A[i][j] = 0.0;
      }
      Y[i] = 0.0;
    }

    for(int i=1; i<nx-1; i++)     // or ia=0; ia<dof; ia++
    {
      x = xg[i];                  // or x = xg[ia+1]
      ia = i-1;
      Y[ia] += f(x); // RHS

      // term  1 * Uxx
      for(int j=0; j<nx; j++)
      {
	if(j==0 || j==nx-1) Y[ia] -= cxx[i][j]*ug[j];
        else A[ia][j-1] += cxx[i][j];
      }

      // term  2 * Ux
      for(int j=0; j<nx; j++)
      {
	if(j==0 || j==nx-1) Y[ia] -= 2.0*cx[i][j]*ug[j];
        else A[ia][j-1] += 2.0*cx[i][j];
      }

      // term  U^2
      // nonlinear section
        A[ia][dof+ia] += 1.0; // just U
  
      // debug
      System.out.println("i="+i+", ia="+ia+", x="+x);
      for(int j=0; j<nxl; j++)
         System.out.println("A["+ia+"]["+j+"]="+A[ia][j]);
      System.out.println("Y["+ia+"]="+Y[ia]);
      System.out.println(" ");
    } // end x

    // debug test
    us[0] = ub(xg[1]); // initial guess, exact solution
    us[1] = ub(xg[2]);
    us[2] = ub(xg[3]);
    SN.simeq(dof, nlin, A, Y, var1, var2, var3, vari1, vari2, us);


    for(int j=0; j<dof; j++) us[j] = 1.0; // initial guess
    SN.simeq(dof, nlin, A, Y, var1, var2, var3, vari1, vari2, us);
    for(int i=1; i<nx-1; i++) ug[i] = us[i-1];

  } // end solve_pde

  void check_soln() // 3 checks, for all interior points
  {
    double err, sumerr, sumsqerr, rmserr, maxerr, avgerr;
    double x, ux, uxx;
    double term, eval_deriv, aterm, aderiv, eterm, ederiv;
    double uAx, uAxx; // numeric derivative, actual solution
    int cnt;

    System.out.println("check solution against PDE");
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;

    for(int i=1; i<nx-1; i++) // inside boundary
    {
      x = xg[i];

      // compute derivatives numerically for this x,y
      // exact ub(xg[i])=uA[i] or computed solution ug[i]

      ux = 0.0;
      uAx = 0.0;
      for(int j=0; j<nx; j++)
      {
        ux += cx[i][j]*ug[j];
        uAx += cx[i][j]*uA[j];
      }
      uxx = 0.0;
      uAxx = 0.0;
      for(int j=0; j<nx; j++)
      {
	uxx += cxx[i][j]*ug[j];
	uAxx += cxx[i][j]*uA[j];
      }

      eval_deriv = 0.0; // based on numeric derivatives using ug[]
      ederiv = 0.0;     // based on numeric derivatives using uA[]
      aderiv = 0.0;     // based on exact derivatives using ub()
      // term   Uxx(x) +
      term = uxx;
      eval_deriv += term;
      eterm = uAxx;
      ederiv += eterm;
      aterm = ubxx(x);
      aderiv += aterm;
      // term   2*Uxx + 
      term = 2.0*ux; 
      eval_deriv += term;
      eterm = 2.0*uAx; 
      ederiv += eterm;
      aterm = 2.0*ubx(x); 
      aderiv += aterm;
      // term  U(x)^2      
      term = ug[i]*ug[i];
      eval_deriv += term;
      eterm = uA[i]*uA[i];
      ederiv += eterm;
      aterm = uA[i]*uA[i];
      aderiv += aterm;

      err = eval_deriv - f(x);
      err = Math.abs(err);
      sumerr += err;
      sumsqerr += err*err;
      maxerr = Math.max(err,maxerr);
      cnt++;  
      // debug
      System.out.println("i="+i+", x="+x+", u="+ub(x));
      System.out.println(" ux="+ux+    ", uxx= "+uxx);
      System.out.println("ubx="+ubx(x)+", ubxx="+ubxx(x));
      System.out.println("uAx="+uAx+   ", uAxx="+uAxx);
      System.out.println("eval="+(eval_deriv)+", f(x)="+f(x));
      System.out.println("err="+err);
      System.out.println("ederiv="+(ederiv)+", f(x)="+f(x));
      System.out.println("eerr="+(ederiv-f(x)));
      System.out.println("aderiv="+(aderiv)+", f(x)="+f(x));
      System.out.println("aerr="+(aderiv-f(x)));
      System.out.println(" ");
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln 

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

} // end class pde_nl21