// simeq_newton2.java  solve nonlinear system of equations
//                     method: newton iteration using Jacobian
//
//  Given problem  A X = Y  where X may have terms x1, x2, x3, x4,
//                          x1*x2, x1*x3, x1*x4, x2*x3, x2*x4, x3*x4
//                 A is 10 by 1 matrix of reals (could be complex)
//                 Y is vector of reals (could be complex)
//                 independent unknowns are x1, x2, x3, x4
//
//  for testing, generate A using pseudo random numbers
//               choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute products
//               compute terms of Y using Y = A X
//
//  Solve by initial guess at values of x1, x2, x3, x4 computing products
//    X_next = X_initial - J_initial^-1 * (A *  X_initial - Y)
//    in general X_next = X_prev - (J_prev^-1 * (A * X_prev - Y))*b
//                                 where 0 < b < 1, often 0.5, for stability
//
//  solved when  abs sum each row A * X_next -Y < epsilon
//
//  It may stall, stop if abs(X_next-X_prev)<epsilon
//  It may diverge, stop, indicate no solution
//                        (or try a different initial guess)
//  It may oscillate, stop, indicate no solution
//                        (or try a different initial guess)
//
//
//  The matrix equation for this case is:
//
//   x1     x2     x3     x4     x1*x2 ... x3*x4 
//
// | A1,1   A1,2   A1,3   A1,4   A1,5 ...  A1,10  | | x1    | |y1 |
// | A2,1   A2,2   A2,3   A2,4   A2,5 ...  A2,10  | | x2    | |y2 |
// | A3,1   A3,2   A3,3   A3,4   A3,5 ...  A3,10  | | x3    | |y3 |
// | A4,1   A4,2   A4,3   A4,4   A4,5 ...  A4,10  |*| x4    |=|y4 |
// | A5,1   A5,2   A5,3   A5,4   A5,5 ...  A5,10  | | x1*x2 | |y5 |
// | A6,1   A6,2   A6,3   A6,4   A6,5 ...  A6,10  | | x1*x3 | |y6 |
// | ...                                          | |       | |   |
// | A10,1  A10,2  A10,3  A10,4  A10,5 ... A10,10 | | x3*x4 | |y10|
//
//  The Jacobian, J is the numerically computed derivatives based on A, X_prev
//  The derivative coefficients may be computed once based on A and the
//  unknown variables. The numeric value plugs in the X_prev values.
//
//  J1,1 = A1,1 + A1,5*x2 + A1,6*x3 + A1,7*x4   deriv Row 1 wrt x1
//  J1,2 = A1,2 + A1,5*x1 + A1,8*x3 + A1,9*x4   deriv Row 1 wrt x2
//  J1,3 = A1,3 + A1,6*x1 + A1,8*x2 + A1,10*x4  deriv Row 1 wrt x3
//  J1,4 = A1,4 + A1,7*x1 + A1,9*x2 + A1,10*x3  deriv Row 1 wrt x4
//  J1,5 = A1,5                                 deriv Row 1 wrt x1, wrt x2
//  J1,6 = A1,6                                 deriv Row 1 wrt x1, wrt x3
//  J1,7 = A1,7                                 deriv Row 1 wrt x1, wrt x4
//  J1,8 = A1,8                                 deriv Row 1 wrt x2, wrt x3
//  J1,9 = A1,9                                 deriv Row 1 wrt x2, wrt x4
//  J1,10 = A1,10                               deriv Row 1 wrt x3, wrt x4
//
//  etc.
//
//  Code or a data structure must be available to know the equation
//  of the entries in the Y vector. Symbolic computation of
//  derivatives is assumed to be available, when needed.
//

public class simeq_newton2
{
  int n = 10; // number of equations
  double A[][] = new double[n][n];
  double Ja[][] = new double[n][n]; // inverted in place
  double Y[] = new double[n];
  double X_prev[] = new double[n];
  double X_temp[] = new double[n];
  double X_tmp2[] = new double[n];
  double X_next[] = new double[n];
  double X_soln[] = new double[n]; // usually not known. test case

  // derivatives hand coded to load Ja

  public simeq_newton2()
  {
    double eps = 1.0e-8;
    double b = 1.00; // stability factor
    double resid;  // residual from last iteration
    double presid; // residual from prior to last iteration
    double err;      // error from test solution

    System.out.println("simeq_newton2 running");
    // build test case
    for(int i=0; i<n; i++) // all subscripts zero based, one less than comments
    {
      for(int j=0; j<n; j++)
      {
	  A[i][j] = Math.random();
      }
    }
    System.out.println("A generated ");
    // choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, x3*x4
    X_soln[0] = 1.1;
    X_soln[1] = 1.2;
    X_soln[2] = 1.4;
    X_soln[3] = 1.5;
    X_soln[4] = X_soln[0]*X_soln[1];
    X_soln[5] = X_soln[0]*X_soln[2];
    X_soln[6] = X_soln[0]*X_soln[3];
    X_soln[7] = X_soln[1]*X_soln[2];
    X_soln[8] = X_soln[1]*X_soln[3];
    X_soln[9] = X_soln[2]*X_soln[3];
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("X_soln["+i+"]="+X_soln[i]);
    }
    System.out.println(" ");

    // set up Y
    for(int i=0; i<n; i++)
    {
      Y[i] = 0.0;
      for(int j=0; j<n; j++)
      {
	Y[i] += A[i][j]*X_soln[j];
      }
    }
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("Y["+i+"]="+Y[i]);
    }
    System.out.println(" ");

    // setup

    // initial condition
    for(int i=0; i<n; i++) X_prev[i] = 1.0;
    // debug print
    for(int i=0; i<n; i++)
    {
      System.out.println("X_prev["+i+"]="+X_prev[i]);
    }
    System.out.println(" ");

    presid = 1.0e12;
    // iterate to find solution
    for(int itr=0; itr<10; itr++)
    {
      // compute residual
      resid = 0.0; 
      for(int i=0; i<n; i++)
      {
        X_temp[i] = 0.0;
        for(int j=0; j<n; j++)
        {
	  X_temp[i] += A[i][j]*X_prev[j];
        }
        X_temp[i] -= Y[i]; // X_temp used later if not solved
        resid += Math.abs(X_temp[i]);
      }
      // debug print
      for(int i=0; i<n; i++)
      {
	System.out.println("residual X_temp["+i+"]="+X_temp[i]);
      }
      System.out.println(" ");
      // for testing, compute error with solution
      err = 0.0;
      for(int i=0; i<n; i++)
      {
	  err += Math.abs(X_prev[i] - X_soln[i]);
      }
      System.out.println("itr "+itr+", prev="+presid+" residual="+resid+
                         ", err="+err);
      // check for convergence
      if(resid<eps) break;

      // compute Jacobian
      for(int i=0; i<n; i++) // all subscripts zero based
      {
        for(int j=0; j<n; j++)
        {
	  Ja[i][j] = A[i][j];
          if(j==0) Ja[i][j] += A[i][4]*X_prev[1] + A[i][5]*X_prev[2] +
		               A[i][6]*X_prev[3];
          if(j==1) Ja[i][j] += A[i][4]*X_prev[0] + A[i][7]*X_prev[2] +
                               A[i][8]*X_prev[3];
          if(j==2) Ja[i][j] += A[i][5]*X_prev[0] + A[i][7]*X_prev[1] +
                               A[i][9]*X_prev[3];
          if(j==3) Ja[i][j] += A[i][6]*X_prev[0] + A[i][8]*X_prev[1] +
                               A[i][9]*X_prev[2];
        }
      }
      System.out.println("Ja computed ");
      // invert Ja
      new invert(Ja);
      System.out.println("Ja inverted ");

      // X_next = X_prev - (J_prev^-1 * (A X_prev - Y))*b
      // X_next = X_prev - (J_prev^-1 * (X_temp))*b
      // X_next = X_prev - (X_tmp2)*b

      for(int i=0; i<n; i++)
      {
	X_tmp2[i] = 0.0;
        for(int j=0; j<n; j++)
        {
	  X_tmp2[i] += Ja[i][j]*X_temp[j]; 
        }
      }
      // debug print
      for(int i=0; i<n; i++)
      {
	System.out.println("change X_tmp2["+i+"]="+X_tmp2[i]);
      }
      System.out.println(" ");

      for(int i=0; i<n; i++)
      {
        X_next[i] = X_prev[i] - b*X_tmp2[i];
      }
      for(int i=0; i<n; i++)
      {
        System.out.println("X_next["+i+"]="+X_next[i]);
      }
      System.out.println(" ");

      // iterate
      presid = resid; 
      for(int i=0; i<n; i++)
      {
	  X_prev[i] = X_next[i];
      }
      System.out.println(" ");

    } // end iteration    
    System.out.println("simeq_newton2 finished");
  }

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