// pde11p_eq.java   small parts for many points
//                  compile    javac -cp . pde11p_eq.java
//                  execute    java  -cp . pde11p_eq > pde11p_eq_java.out
//
// solution         U(x) = x^4 + 2*x^3 + 3*x^2 + 4*x + 5
// 1st derivative   Ux(x) = 4*x^3 + 6*x^2 + 6*x + 4
// 2nd derivative   Uxx(x) = 12*x^2  + 12*x + 6
// ode since only one variable  Uxx + 2*U = f(x)
// f(x) = Uxx(x) + 2*U(x)
// f(x) = 2*x^4 + 4*x^3 + 18*x^2 + 20*x + 16 
// xmin = -1.0, xmax = 1.0, boundary conditions  f(-1.0), f(1.0)

import java.text.*;
import java.io.*;

public class pde11p_eq
{

  int nx = 11;        // all of grid includes boundary 0 .. nx-1
  int neqn = nx-2;    // number of simultaneous equations to solve offset S(i)
  int npmax = 7;      // points for computing numerical derivative by parts

  boolean debug = true;
  double xg[] = new double[nx];
  double cxx[] = new double[nx];
  double xmin = -1.0;
  double xmax = +1.0;
  double hx;
  double error;         // sum of absolute exact-computed
  double avgerror;      // average error
  double rmserror;      // root mean square error L2 norm       
  double maxerror;      // maximum cell error Linf norm         

  DecimalFormat fe = new DecimalFormat("0.00E00");
  DecimalFormat f6 = new DecimalFormat("0.00000");

  public pde11p_eq()
  {
    System.out.println("pde11p_eq.java running");
    System.out.println("uniform grid on line -1 to 1");
    System.out.println("PDE to solve Uxx + 2*Ux + 3*U = f(x)");
    System.out.println("Uxx(x) + 2*U(x) - 3 = f(x)");
    System.out.println("f(x) = 2*x^4 + 4*x^3 + 18*x^2 + 20*x + 13");
    System.out.println("known function, for testing method");
    System.out.println("U(x) := x^4 + 2*x^3 + 3*x^2 + 4*x + 5");

    System.out.println(" ");

    hx = (xmax-xmin)/(double)(nx-1);
    System.out.println("xmin="+xmin+", xmax="+xmax+", nx="+nx+", hx="+hx);

    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      System.out.println("xg("+i+")="+xg[i]+", U(xg[i])="+U(xg[i]));
      System.out.println("                           f(xg[i])="+F(xg[i]));
    }

    System.out.println("xmin="+f6.format(xg[0])+
		       ", xmax="+f6.format(xg[nx-1])+
                       ", nx="+nx+", hx="+f6.format(hx)); 
    System.out.println("U(xmin)(xmax)    ="+f6.format(U(xmin))+"  "+
                                      f6.format(U(xmax)));
    System.out.println("Uxx(xmin)(xmax)  ="+f6.format(Uxx(xmin))+"  "+
                                      f6.format(Uxx(xmax)));
    System.out.println("f(xmin)(xmax)    ="+f6.format(F(xmin))+"  "+
                                      f6.format(F(xmax)));
    System.out.println(" ");
    System.out.println("check normal PDE method ");
    solve_pde();
    System.out.println(" ");
    System.out.println("check 3 part PDE method ");
    solve_pdep(); // 3 parts
    System.out.println(" ");
    System.out.println("pde11p_eq.java finished");
  }

  void  solve_pdep() // in 3 parts for each order derivative
  {
    int np = 7;      // local, should test 5,7,9
    int iv = 0;      // derivative point and RHS f(xg[S(iv)]

    double A[][] = new double[neqn][neqn+1]; // last column is RHS 
    double Us[] = new double[neqn];
    double Uxxval, eqnxxval;
    double val, sum, err;
    int row;
    int cs = neqn;

    System.out.println(" ");
    System.out.println(" ");
    System.out.println("solve_pdep using 3 part method nx="+nx+", neqn="+neqn+
                       ", np="+np);

    // zero A
    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<=neqn; j++) A[i][j] = 0.0;
    }

    for(int i=1; i<nx-1; i++) // S(i) 0..neqn-1  0..nx-3
    {
      row = S(i); // not including boundary, known
      System.out.println("row = "+row);

      // for each term in PDE

      double coefdxx = 1.0;
      if(i<(np+1)/2) // first section
      {
        System.out.println("i="+i+", (np+1)/2="+((np+1)/2));
        new rderiv(2,np,i,hx,cxx);
        if(debug) System.out.println("1 rderiv(2,np,"+i+",hx,cxx)");
        for(int ip=0; ip<np; ip++)
        {
	  iv=0;
	  System.out.println("xx1 i="+i+", ip="+ip+", iv="+iv);
          val = coefdxx*cxx[ip];
          if(ip==0) // boundary
	  {
            A[row][cs] -= val*U(xg[0]);
	  }
          else
	  {
	    A[row][S(ip)] += val;
            System.out.println("xx1 A["+row+"]["+S(ip)+"] added cxx["+ip+"]");
	  }
        } // end ip

        // for U(x)                                                   
        double coefu = 2.0;
        A[row][S(i)] += coefu;

        A[row][cs] += F(xg[i]);

        if(debug) // debug
        {
          sum = 0.0;
          for(int j=0; j<neqn; j++)
	      sum += A[row][j]*U(xg[j+1]);
	  err = sum - A[row][cs];
          System.out.println("sum="+sum+", row error="+err);
          System.out.println(" ");
	}
      } // end first

      if(i>=(np+1)/2 && i<nx-(np+1)/2) // second section no boundary
      {
	new rderiv(2,np,(np-1)/2,hx,cxx);
        if(debug) System.out.println("xx2 rderiv(2,np,"+((np-1)/2)+",hx,cxx)");
        iv = i-(np+1)/2;
        for(int ip=0; ip<np; ip++) // no boundary
        {
	  System.out.println("xx2 i="+i+", ip="+ip+", iv="+iv);
          val = coefdxx*cxx[ip];
	  A[row][ip+iv] += val;
          System.out.println("xx2 A["+row+"]["+(ip+iv)+"] added cx["+ip+"]");
        }
        // for U(x)                                                   
        double coefu = 2.0;
        A[row][S(i)] += coefu;

        A[row][cs] += F(xg[i]);

        if(debug) // debug
        {
          sum = 0.0;
          for(int j=0; j<neqn; j++)
	      sum += A[row][j]*U(xg[j+1]);
	  err = sum - A[row][cs];
          System.out.println("sum="+sum+", row error="+err);
          System.out.println(" ");
	}
      } // end second

      if(i>=nx-(np+1)/2) // third section  upper boundary
      {
	new rderiv(2,np,np-(nx-i),hx,cxx);
        if(debug) System.out.println("xx3 rderiv(2,np,"+(np-(nx-i))+",hx,cxx)");
        iv = i-(np-(nx-i))-1;
        for(int ip=0; ip<np; ip++)
        { 
	  System.out.println("xx3 i="+i+", ip="+ip+", iv="+iv);
          val = coefdxx*cxx[ip];
          if(ip==np-1)
	  {
	    A[row][cs] -= val*U(xg[nx-1]); // upper boundary
	  }
          else
	  {
	    A[row][ip+iv] += val;
            System.out.println("xx3 A["+row+"]["+(ip+iv)+"] added cxx["+ip+"]");
	  }
	} // end ip

        // for U(x)                    
        double coefu = 2.0;
        A[row][S(i)] += coefu;

        // for f(x) RHS constant    
        A[row][cs] += F(xg[i]);

        if(debug) // debug
        {
          sum = 0.0;
          for(int j=0; j<neqn; j++)
	      sum += A[row][j]*U(xg[j+1]);
	  err = sum - A[row][cs];
          System.out.println("sum="+sum+", row error="+err);
          System.out.println(" ");
	}
      } // end third
    } // end i
    // if(debug) prt_mat(neqn, A);
    check_rows(neqn, A);

    // solve for U
    new simeq(neqn, A, Us);

    // check for error
    error = 0.0;
    avgerror = 0.0;      // average error
    rmserror = 0.0;      // root mean square error L2 norm       
    maxerror = 0.0;      // maximum cell error Linf norm
    for(int i=0; i<neqn; i++)
    {
      error = Math.abs(Us[i]-U(xg[i+1])); // not boundary
      avgerror = avgerror + error;
      rmserror = rmserror + error*error;
      maxerror = Math.max(maxerror, error);
      System.out.println("i="+i+", Us[i]="+Us[i]+", U(xg[i+1])="+U(xg[i+1]));
    }
    avgerror = avgerror/(double)neqn;
    rmserror = Math.sqrt(rmserror/(double)neqn);
    System.out.println("avg err="+avgerror);
    System.out.println("rms err="+rmserror);
    System.out.println("max err="+maxerror);

  } // end solve_pdep

  void  solve_pde() // all in one standard method
  {
    double A[][] = new double[neqn][neqn+1]; // last column is RHS 
    double Us[] = new double[neqn];
    double Uxval, Uxxval;
    double val, sum, err;
    int row;
    int cs = neqn;

    System.out.println("solve_pde using  nx="+nx+", neqn="+neqn);

    // zero A
    for(int i=1; i<nx-1; i++)
    {
	for(int j=1; j<nx; j++) A[S(i)][S(j)] = 0.0;
    }

    for(int i=1; i<nx-1; i++)
    {
      row = S(i); // not including boundary, known
      System.out.println("row = "+row);
      // for each term in PDE

      // for d^2u/dx^2                                                    
      double coefdxx = 1.0;
      new rderiv(2,nx,i,hx,cxx);
      for(int j=0; j<nx; j++)
      {
        val = coefdxx*cxx[j];
        if(j==0 || j==nx-1)
	{
	  A[row][cs] -= val*U(xg[j]);
	}
        else
	  A[row][S(j)] += val;
      }

      // for U(x)                                                   
      double coefu = 2.0;
      A[row][S(i)] += coefu;

      // for f(x) RHS constant
      A[row][cs] += F(xg[i]);
      if(row<3) // debug
      {
        sum = 0.0;
        for(int j=0; j<neqn; j++)
	    sum += A[row][j]*U(xg[j+1]);
	err = sum - A[row][cs];
        // System.out.println("sum="+sum+", row error="+err);
        // System.out.println(" ");
      }

    } // end i
    // if(debug) prt_mat(neqn, A);
    check_rows(neqn, A);

    // solve for U
    new simeq(neqn, A, Us);

    // check for error
    error = 0.0;
    avgerror = 0.0;      // average error
    rmserror = 0.0;      // root mean square error L2 norm       
    maxerror = 0.0;      // maximum cell error Linf norm
    for(int i=0; i<neqn; i++)
    {
      error = Math.abs(Us[i]-U(xg[i+1])); // not boundary
      avgerror = avgerror + error;
      rmserror = rmserror + error*error;
      maxerror = Math.max(maxerror, error);
      System.out.println("i="+i+", Us[i]="+Us[i]+", U(xg[i+1])="+U(xg[i+1]));
    }
    avgerror = avgerror/(double)neqn;
    rmserror = Math.sqrt(rmserror/(double)neqn);
    System.out.println("avg err="+avgerror);
    System.out.println("rms err="+rmserror);
    System.out.println("max err="+maxerror);

  } // end solve_pde

  double U(double x) // x^4 + 2*x^3 + 3*x^2 + 4*x + 5
  {
    return x*x*x*x + 2.0*x*x*x + 3.0*x*x + 4.0*x + 5.0; 
  } // end U

  double Uxx(double x) //  12*x^2  + 12*x + 6
  {
    return  12.0*x*x + 12.0*x + 6.0; 
  } // end Uxx

  double F(double x)
  {
    return  2.0*x*x*x*x + 4.0*x*x*x + 18.0*x*x + 20.0*x + 16.0 ;
  } // end F

  int S(int i)
  {
    return i-1;
  } // end S

    void prt_mat(int ne, double A[][])
  {
    for(int i=0; i<ne; i++)
      for(int j=0; j<=ne; j++)
        System.out.println("A["+i+"]["+j+"]="+A[i][j]);
  } // end prt_mat

    void check_rows(int ne, double A[][]) // A[ne][ne+1]
  {
    double val = 0.0;
    System.out.println("check rows for equal zero");
    for(int i=0; i<ne; i++)
    {
      val = 0.0;
      for(int j=0; j<ne; j++)
      {
	val += A[i][j]*U(xg[j+1]);
      }
      val = val - A[i][ne]; // should be near zero
      System.out.println("row="+i+", err="+val+", RHS="+A[i][ne]);
    }
    System.out.println(" ");
  } // end check_rows

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