// pde_i1.java   initial value PDE  very simple test cases
// find values of  U(x,y) given only U[0,0] 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(0,0) =1  given
// see http://www.cs.umbc.edu/~squire/cs455_l27b.shtml  for development

public class pde_i1
{
  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)
  {
    return (Ux(x,y) + Ux(x+h/2.0,y) + Ux(x+h,y))/3.0;
  }

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

  public pde_i1()
  {
    double U[][] = new double[5][5];
    double hx = 0.1;
    double hy = 0.15;
    double x, y;
    double xp, yp, err, tmp;
    double dx0, dx1, dx2, dx3, dy0, dy1, dy2, dy3;

    System.out.println("pde_i1.java running");
    System.out.println(" ");
    System.out.println("some preliminary derivative checks");
    dx0 = (-3.0*U(0.0,0.0) + 4.0*U(hx,0.0) - 1.0*U(2.0*hx,0.0))/(2.0*hx);
    err = dx0 - Ux(0.0,0.0);
    System.out.println("Ux0="+Ux(0.0,0.0)+", dx0="+dx0+", err="+err);
    dx1 = (-1.0*U(0.0,0.0) + 1.0*U(2.0*hx,0.0))/(2.0*hx);
    err = dx1 - Ux(hx,0.0);
    System.out.println("Ux1="+Ux(hx,0.0)+", dx1="+dx1+", err="+err);
    dx2 = (1.0*U(0.0,0.0) - 4.0*U(hx,0.0) + 3.0*U(2.0*hx,0.0))/(2.0*hx);
    err = dx2 - Ux(2.0*hx,0.0);
    System.out.println("Ux2="+Ux(2.0*hx,0.0)+", dx2="+dx2+", err="+err);
    System.out.println(" ");
    dy0 = (-3.0*U(0.0,0.0) + 4.0*U(0.0,hy) - 1.0*U(0.0,2.0*hy))/(2.0*hy);
    err = dy0- Uy(0.0,0.0);
    System.out.println("Uy0="+Uy(0.0,0.0)+", dy0="+dy0+", err="+err);
    dy1 = (-1.0*U(0.0,0.0) + 1.0*U(0.0,2.0*hy))/(2.0*hy);
    err = dy1- Uy(0.0,hy);
    System.out.println("Uy1="+Uy(0.0,hy)+", dy1="+dy1+", err="+err);
    dy2 = (1.0*U(0.0,0.0) - 4.0*U(0.0,hy) + 3.0*U(0.0,2.0*hy))/(2.0*hy);
    err = dy2- Uy(0.0,2.0*hy);
    System.out.println("Uy2="+Uy(0.0,2.0*hy)+", dy2="+dy2+", err="+err);
    System.out.println(" ");

    dx1 = (-3.0*U(hx,hy) + 4.0*U(2.0*hx,hy) - 1.0*U(3.0*hx,hy))/(2.0*hx);
    err = dx1 - Ux(hx,hy);
    System.out.println("Ux1="+Ux(hx,hy)+", dx1="+dx1+", err="+err);
    dx2 = (-1.0*U(hx,hy) + 1.0*U(3.0*hx,hy))/(2.0*hx);
    err = dx2 - Ux(2.0*hx,hy);
    System.out.println("Ux2="+Ux(2.0*hx,0.0)+", dx2="+dx2+", err="+err);
    dx3 = (1.0*U(hx,hy) - 4.0*U(2.0*hx,hy) + 3.0*U(3.0*hx,hy))/(2.0*hx);
    err = dx3 - Ux(3.0*hx,hy);
    System.out.println("Ux3="+Ux(3.0*hx,hy)+", dx3="+dx3+", err="+err);
    System.out.println(" ");
    dy1 = (-3.0*U(hx,hy) + 4.0*U(hx,2.0*hy) - 1.0*U(hx,3.0*hy))/(2.0*hy);
    err = dy1- Uy(hx,hy);
    System.out.println("Uy1="+Uy(hx,hy)+", dy1="+dy1+", err="+err);
    dy2 = (-1.0*U(hx,hy) + 1.0*U(hx,3.0*hy))/(2.0*hy);
    err = dy2- Uy(hx,2.0*hy);
    System.out.println("Uy2="+Uy(hx,2.0*hy)+", dy2="+dy2+", err="+err);
    dy3 = (1.0*U(hx,hy) - 4.0*U(hx,2.0*hy) + 3.0*U(hx,3.0*hy))/(2.0*hy);
    err = dy3- Uy(hx,3.0*hy);
    System.out.println("Uy3="+Uy(hx,3.0*hy)+", dy3="+dy3+", err="+err);
    System.out.println(" ");


    U[0][0] = 1.0;
    System.out.println("hx="+hx+", hy="+hy+", U[0][0]="+U[0][0]);
    System.out.println(" ");
    x = 0.0;
    y = 0.0;

    System.out.println("using simple U(x+h,y) = U(x,y)+hx*Ux(x,y)");
    U[1][0] = U[0][0] + hx*Ux(x,y);
    err = U[1][0] - U(x+hx,y);
    System.out.println("U[1][0] = U(x+hx,y)="+U[1][0]+", err="+err+" bad");
    U[0][1] = U[0][0] + hy*Uy(x,y);
    err = U[0][1] - U(x,y+hy);
    System.out.println("U[0][1] = U(x,y+hy)="+U[0][1]+", err="+err+" bad");

    System.out.println(" ");
    System.out.println("using simple U(x+h,y) = U(x,y)+hx*Ux(x+hx,y)");
    U[1][0] = U[0][0] + hx*Ux(x+hx,y);
    err = U[1][0] - U(x+hx,y);
    System.out.println("U[1][0] = U(x+hx,y)="+U[1][0]+", err="+err+" bad");
    U[0][1] = U[0][0] + hy*Uy(x,y+hy);
    err = U[0][1] - U(x,y+hy);
    System.out.println("U[0][1] = U(x,y+hy)="+U[0][1]+", err="+err+" bad");

    System.out.println(" ");
    System.out.println("using U(x+h,y) = U(x,y)+better derivative");
    System.out.println("exact U[1][0]="+U(x+hx,y)+", U[2][0]="+U(x+2.0*hx,y)+
                       ", U[3][0]="+U(x+3.0*hx,y));
    U[1][0] = U[0][0] + hx*dx(x,y,hx);
    err = U[1][0] - U(x+hx,y);
    System.out.println("U[1][0] = U(x+hx,y)="+U[1][0]+", err="+err);
    x = x +hx;
    y = 0.0;
    U[2][0] = U[1][0] + hx*dx(x,y,hx);
    err = U[2][0] - U(x+hx,y);
    System.out.println("U[2][0] = U(x+2*hx,y)="+U[2][0]+", err="+err);
    x = x +hx;
    U[3][0] = U[2][0] + hx*dx(x,y,hx);
    err = U[3][0] - U(x+hx,y);
    System.out.println("U[3][0] = U(x+3*hx,y)="+U[3][0]+", err="+err);

    System.out.println(" ");
    System.out.println("exact U[0][1]="+U(x,y+hy)+", U[0][2]="+U(x,y+2.0*hy)+
                       ", U[0][3]="+U(x,y+3.0*hy));
    x = 0.0;
    y = 0.0;
    U[0][1] = U[0][0] + hy*dy(x,y,hy);
    err = U[0][1] - U(x,y+hy);
    System.out.println("U[0][1] = U(x,hy)="+U[0][1]+", err="+err);
    y = y+hy;
    U[0][2] = U[0][1] + hy*dy(x,y,hy);
    err = U[0][2] - U(x,y+hy);
    System.out.println("U[0][2] = U(x,y+2.0*hy)="+U[0][2]+", err="+err);
    y = y+hy;
    U[0][3] = U[0][2] + hy*dy(x,y,hy);
    err = U[0][3] - U(x,y+hy);
    System.out.println("U[0][3] = U(x,y+3.0*hy)="+U[0][3]+", err="+err);

    System.out.println(" ");
    System.out.println("exact U[1][1]="+U(hx,hy)+", U[2][2]="+U(2.0*hx,2.0*hy));
    x = hx;  // initial point [1][0]
    y = 0.0;
    U[1][1] = U[1][0] + hy*dy(x,y,hy); // moving up
    tmp = U[1][1];
    err = U[1][1] - U(hx,hy);
    System.out.println("U[1][1] = U(hx,hy)="+U[1][1]+", yerr="+err);
    x = 0.0;
    y = hy;
    U[1][1] = U[0][1] + hx*dx(x,y,hx); // moving right
    err = U[1][1] - U(hx,hy);
    System.out.println("U[1][1] = U(hx,hy)="+U[1][1]+", xerr="+err);
    tmp = (tmp + U[1][1])/2.0;
    err = tmp - U(hx,hy);
    System.out.println("U[1][1] = U(hx,hy)="+tmp+", avgerr="+err);

    System.out.println(" ");
    U[1][2] = U(hx,2.0*hy);
    U[2][1] = U(2.0*hx,hy);
    x = hx;
    y = 2.0*hy;
    U[2][2] = U[1][2] + hx*dx(x,y,hx); // moving right
    tmp = U[2][2];
    err = U[2][2] - U(2.0*hx,2.0*hy);
    System.out.println("U[2][2] = U(2.0*hx,2.0*hy)="+U[2][2]+", xerr="+err);
    x = 2.0*hx;
    y = hy;
    U[2][2] = U[2][1] + hy*dy(x,y,hy); // moving up
    err = U[2][2] - U(2.0*hx,2.0*hy);
    System.out.println("U[2][2] = U(2.0*hx,2.0*hy)="+U[2][2]+", yerr="+err);
    tmp = (tmp + U[2][2])/2.0;
    err = tmp - U(2.0*hx,2.0*hy);
    System.out.println("U[2][2] = U(2.0*hx,2.0*hy)="+tmp+", avgerr="+err);
    System.out.println(" ");

    System.out.println("pde_i1.java finished");
  }

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