// test_pde_nl.java       compile    javac -cp . test_pde_nl.java
//                        execute    java  -cp . test_pde_nl
// see test_nav2_newton3.java
// needs simeq_newton3.java
// U(x) + V(x) + Ux(x) + Vx(x) + U(x)*Ux(x) + V(x)*Vx(x) = F(x)
// U(x) = 2*x*x = x1
// V(x) = 3*(x+1)*(x+1)x = x2
// Ux(x) = 4*x  = x3
// Vx(x) = 6*x+6  = x4
public class test_pde_nl
{
  int neqn = 4;
  int nlin = 2;
  int ntot = neqn+nlin;
  double A[][] = new double[neqn][ntot];
  //
  //coef U(x)   V(x)   Ux(x)  Vx(x)  U*Ux   V*Vx      unk           Y
  //   | A1,1   A1,2   A1,3   A1,4   A1,5   A1,6  | | U(x)       | |Y1 |
  //   | A2,1   A2,2   A2,3   A2,4   A2,5   A2,6  | | V(x)       | |Y2 |
  //   | A3,1   A3,2   A3,3   A3,4   A3,5   A3,6  |*| Ux(x)      |=|Y3 |
  //   | A4,1   A4,2   A4,3   A4,4   A4,5   A4,6  | | Vx(x)      | |Y4 |
  //                                              | | U(x)*Ux(x) |
  //                                              | | V(v)*Vx(x) |
  //
  double Y[] = new double[neqn];
  double X[] = new double[ntot]; // solution A X = Y
  int var1[] = new int[ntot]; // neqn variables then nlin
  int var2[] = new int[ntot]; // neqn -1 then nlin
  int var3[] = new int[ntot]; // all -1
  int maxitr = 2; // maximum iterations of newton method
  double b = 0.5; // stability factor
  double eps = 1.0e-4;
  Boolean debug = true;
    
  public test_pde_nl() // constructor
  {
    double x = 1.0; // location of solution
    double err = 0.0;
    System.out.println("test_pde_nl.java running");
    System.out.println("neqn="+neqn+", nlin="+nlin+", ntot="+ntot);
    System.out.println("maxitr="+maxitr+", b="+b+", eps="+eps);
    System.out.println("U(x)+V(x)+Ux(x)+Vx(x)+U(x)*Ux(x)+V(x)*Vx(x)=F(x)");
  
    System.out.println("Making newton3 var1, var2, var3");
    for(int i=0; i<ntot; i++)
    {
       var3[i] = -1; // no v^3
    }
    for(int i=0; i<neqn; i++)
    {
      var1[i] = i;  // just coef
      var2[i] = -1; // no V^2
    }
    // for(int i=neqn; i<ntot; i++)
    {
      var1[4] = 0; // product 
      var2[4] = 2;
      var1[5] = 1; // product 
      var2[5] = 3; // var3 = -1 above
    }

    System.out.println("   var1, var2, var3 ");
    for(int i=0; i<ntot; i++)
    {
	System.out.println("["+i+"] "+var1[i]+",   "+var2[i]+",   "+var3[i]);
    }
    System.out.println(" ");

    init_matrix();

    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<ntot; j++)
      {
	System.out.println("A["+i+"]["+j+"] ="+A[i][j]);
      } // end j
      System.out.println("Y["+i+"] ="+Y[i]);
      System.out.println(" ");
    }
    System.out.println(" ");

    x = 1.0; // initial guess
    // for(int j=0; j<ntot; j++)
    // {
    X[0] = U(x);
    X[1] = V(x);
    X[2] = Ux(x);
    X[3] = Vx(x);
    X[4] = U(x)*Ux(x);
    X[5] = V(x)*Vx(x);      
    // end j
    System.out.println("initial X[0]="+X[0]+"  X[1]="+X[1]);
    System.out.println("initial X[2]="+X[2]+"  X[3]="+X[3]);
    System.out.println("initial X[4]="+X[4]+"  X[5]="+X[5]);
    System.out.println(" ");
    err = check_eqn(neqn, ntot, A, Y, X); // A*X=Y
    System.out.println("check_eqn err= "+err);
    System.out.println(" ");
    
    new simeq_newton3(neqn, nlin, A, Y, var1, var2, var3, X,
		      maxitr, b, eps, 0);
    x = 1.0;
    System.out.println("want X[0]="+U(x)+"  X[1]="+V(x));
    System.out.println("want X[2]="+Ux(x)+"  X[3]="+Vx(x));
    System.out.println("want X[4]="+(U(x)*Ux(x))+"  X[5]="+(V(x)*Vx(x)));
    System.out.println("solution X[0] X[1] X[2] X[3] X[4] X[5] ");
    System.out.println(X[0]+"  "+X[1]+"  "+X[2]+"  "+X[3]+"  "+X[4]+"  "+X[5]);
    System.out.println(" ");

    // check solution

    System.out.println("test_pde_nl.java finished");
  } // end test_pde_nl

  void init_matrix()
  {
    double x = 1.0;
    double xr = 1.0;
    for(int i=0; i<neqn; i++)
    {
      // for(int j=0; j<ntot; j++)
      A[i][0] = U(x)/U(xr); 
      A[i][1] = V(x)/V(xr); 
      A[i][2] = Ux(x)/Ux(xr); 
      A[i][3] = Vx(x)/Vx(xr);
      A[i][4] = (U(x)*Ux(x))/(U(xr)*Ux(xr)); 
      A[i][5] = (V(x)*Vx(x))/(V(xr)*Vx(xr)); 

      // end j
      // Y[i] = F(xr);
      Y[i] = Fc(x, A[i][0],  A[i][1],  A[i][2],  A[i][3],  A[i][4],  A[i][5]);
      xr = xr+0.1; // do not allow det(A) == 0
    } // end i
  } // end init_matrix

  // F(x) = U(x) + v(x) + Ux(x) + Vx(x) + U(x)*Ux(x) + V(x)*Vx(x)
  double F(double x)
  {
    return U(x) + V(x) + Ux(x) + Vx(x) +
	   U(x)*Ux(x) + V(x)*Vx(x);
  } // end F

  // Fc(x) = U(x) + V(x) + Ux(x) + Vx(x) + U(x)*Ux(x) + V(x)*Vx(x)
  double Fc(double x,  double c1, double c2, double c3, double c4,
	   double c5, double c6)
  {
    return c1*U(x) + c2*V(x) + c3*Ux(x) + c4*Vx(x) +
	   c5*U(x)*Ux(x) + c6*V(x)*Vx(x);
  } // end Fc

  double U(double x)   // U(x) = 2*x*x
  {
    return 2.0*x*x;
  } // end U

  double V(double x)   // V(x) = 3*x*x
  {
    return 3.0*(x+1.0)*(x+1.0);
  } // end V

  double Ux(double x) // Ux(x) = 4*x 
  {
    return 4.0*x;
  } // end Ux

  double Vx(double x) // Vx(x) = 6*x
  {
    return 6.0*x+6.0;
  } // end Vx

  double check_eqn(int n, int nt, double A[][], double Y[], double X[])
  {
    double err=0.0;
    double sum = 0.0;
    for(int i=0; i<n; i++)
    {
      sum = 0.0;
      for(int j=0; j<nt; j++)
      {
	sum = sum + A[i][j]*X[j];
      } // end j
      err = Math.max(err,Math.abs(sum-Y[i]));
    } // end i
    return err;
  } // end check_eqn
    
  public static void main (String[] args)
  {
    new test_pde_nl();
  }
} // end test_pde_nl.java