// pde_i2.java   PDE initial condition example
// find values of  U(x,y) given only U[x0,y0] and
//                        given derivative functions for Ux(x,y), Uy(x,y)
// test case unknown, for testing, U(x,y) = x^2 + 2*x*y + 3*y*y + 1
// Ux(x,y) = 2*x + 2*y     Uy(x,y) = 2*x + 6*y  U(x0,y0) =1  given
// x0 = 0, y0 = 0, xmax = 1, ymax =1, nx = 5, ny = 6         given
//
// see http://www.cs.umbc.edu/~squire/cs455_l27b.shtml  for development
// see pde_i1.java for preliminary tests


public class pde_i2
{
  double U(double x, double y) // only used for error test
  {
    return x*x + 2.0*x*y + 3.0*y*y + 1.0;
  }

  double Ux(double x, double y)
  {
    return 2.0*x + 2.0*y;
  }

  double Uy(double x, double y)
  {
    return 2.0*x + 6.0*y;
  }

  double dx(double x, double y, double h) // helper function
  {
    return (Ux(x,y) + Ux(x+h/2.0,y) + Ux(x+h,y))/3.0;
  }

  double dy(double x, double y, double h) // helper function
  {
    return (Uy(x,y) + Uy(x,y+h/2.0) + Uy(x,y+h))/3.0;
  }

  public pde_i2()
  {
    int nx = 5;
    int ny = 6;
    double x0 = 0.0;
    double y0 = 0.0;
    double xmax = 1.0;
    double ymax = 1.0;
    double initial_value = 1.0;

    double U[][] = new double[nx][ny];
    double hx = (xmax-x0)/(double)(nx-1);
    double hy = (ymax-y0)/(double)(ny-1);
    double x, y;
    double err=0.0, sumerr=0.0, maxerr=0.0;
    int i, j;

    System.out.println("pde_i2.java running nx="+nx+", ny="+ny);
    System.out.println("hx="+hx+", hy="+hy);
    System.out.println("xmax="+xmax+", ymax="+ymax);
    U[0][0] = initial_value;
    x = x0;
    y = y0;
    System.out.println("x0="+x0+", y0="+y0+", U[0][0]="+U[0][0]);
    System.out.println(" ");

    for(j=0; j<ny; j++)
    {
      if(j>0) // first already initialized
      {
        x = x0;
        i = 0;
        U[i][j] = U[i][j-1] + hy*dy(x,y-hy,hy);
        err = Math.abs(U[i][j]-U(x,y));
        sumerr = sumerr + err;
        maxerr = Math.max(maxerr, err);
      }
      for(i=0; i<nx-1; i++)
      {
        U[i+1][j] = U[i][j] + hx*dx(x,y,hx);
        err = Math.abs(U[i+1][j]-U(x+hx,y));
        sumerr = sumerr + err;
        maxerr = Math.max(maxerr, err);
        x = x + hx;
      } // end i
      y = y + hy;
    } // end j

    for(j=0; j<ny; j++)
    {
      for(i=0; i<nx; i++)
      {
	System.out.println("U["+i+"]["+j+"]="+U[i][j]);
      }
    }
    System.out.println("maxerr="+maxerr+", avgerr="+(sumerr/(double)(nx*ny)));
    System.out.println("pde_i2.java finished");
  }

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