// pde3.java  just checking second order approximation         
//      u(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  
//       du/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  
//   d^2u/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  
//       du/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  
//   d^2u/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  
//                                                          
// solve d^2u/dx^2 + d^2u/dy^2 = 12x + 10y = c(x,y)         
// by second order method on square -1,-1 to 1,1            

// second order numerical difference 
// d^2u/dx^2 + d^2u/dy^2 = 
// 1/h^2(u(x-h,y)-2u(x,y)+u(x+h,y))+ 
// 1/h^2(u(x,y-h)-2u(x,y)+u(x,y+h)) + O(h^2)
//  = c(x,y) 
// solve for new u(x,y) = 
// (u(x-h,y)+u(x+h,y)+u(x,y-h)+u(x,y+h)-c(x,y)*h^2)/4

public class pde3
{
  final int nx = 11;
  final int ny = 11;
  double us[][] = new double[nx][ny]; // solution being iterated        
  double ut[][] = new double[nx][ny]; // temporary for partial results 
                     // nx-2 by ny-2 internal, non boundary cells 
  final double xmin = -1.0;
  final double ymin = -1.0;
  final double xmax =  1.0;
  final double ymax =  1.0;
  double h;   // x  and  y step = (xmax-xmin)/(n-1) = (ymax-ymin) 
  double diff=1.0e100;  // sum abs difference form last iteration 
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  public pde3()
  {
    int iter=0;
  
    System.out.println("pde3.java running");
    h = (xmax-xmin)/(nx-1); // step size in x 
    //= (ymax-ymin)/(ny-1);    step size in y 
    System.out.println("xmin="+xmin+ "xmax="+xmax+ "nx="+nx+ "h="+h);
    System.out.println("ymin="+ymin+ "ymax="+ymax+ "ny="+ny+ "h="+h);
  
    printexact();
    init_boundary();
    System.out.println(" initial condition");
    printus();
  
    while(diff>0.0001 && iter<200)
    {
      iterate();
      check();
      avgerror = error/((nx-2)*(ny-2));
      System.out.println("iter="+iter+ "  diff="+diff+ 
                         "  error="+error+ "  avg="+avgerror+ "  max="+maxerror);
      iter++;
      if(iter%50==0) { printus(); printdiff(); }
    }
  } // end pde3

  double u(double x, double y) // solution for boundary and check 
  {
    return x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0;
  }

  double c(double x, double y) // from solve 
  {
    return 12.0*x + 10.0*y;
  }

  double u00(int i, int j) // new iterated value 
  {
    return ( us[i-1][j] + us[i+1][j] + us[i][j-1] + us[i][j+1] -
             c(xmin+i*h,ymin+j*h)*h*h )/4.0;
  }

  void init_boundary()
  {
    int i, j;
  
    // set boundary on four sides 
    for(i=0; i<nx; i++)
    {
      us[i][0]    = u(xmin+i*h, ymin);
      us[i][ny-1] = u(xmin+i*h, ymax);
    }
    for(j=0; j<ny; j++)
    {
      us[0][j]    = u(xmin, ymin+j*h);
      us[nx-1][j] = u(xmax, ymin+j*h);
    }
    // zero internal cells 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        us[i][j] = 0.0;
  } // end init_boundary 

  void iterate()
  {
    int i, j;
  
    diff = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        ut[i][j] = u00(i,j); 
        diff = diff + Math.abs(ut[i][j]-us[i][j]);
      }
    }
    // copy temp to solution 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        us[i][j] = ut[i][j];
  } // end iterate

  void check() // typically not known, instrumentation 
  {
    int i, j;
    double uexact, err;
  
    error = 0.0;
    maxerror = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        uexact = u(xmin+i*h,ymin+j*h);
        err = Math.abs(us[i][j]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  } // end check

  void printus()
  {
    int i, j;
  
    System.out.println("\n computed solution");
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        System.out.println("us["+i+"]["+j+"]= "+us[i][j]);
      }
    }
  } // end printus

  void printexact()
  {
    int i, j;
  
    System.out.println("\n exact solution");
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        System.out.println("u["+i+"]["+j+"]= "+u(xmin+i*h,ymin+j*h));
      }
    }
  } // end printexact

  void printdiff()
  {
    int i, j;
  
    System.out.println("\n difference exact-computed solution");
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        System.out.println("diff["+i+"]["+j+"]= "+(u(xmin+i*h,ymin+j*h)-us[i][j]));
      }
    }
  } // end printdiff
  
  public static void main (String[] args)
  {
    new pde3();
  } // end main
} // end class pde3 of pde3.java