// rk1_xy_acc.java   1st order ODE  dyf(x)/dx = yp(x) = 4*x^3  given
//               check solution yf(x) =  x^4
// the initial conditions x0=0 y=0, xn=2, h=0.1 solve for y(y) = y^4
// Runge Kutta 4th order, third order, second order, first order

public class rk1_xy_acc
{
  public rk1_xy_acc()
  {
    double x0 = 1.0;   // starting x
    double xn = 3.0;   // ending   x uses parts 
    double y0  = 1.0;  // y at x0  or starting x
    double h  = 0.1;   // step size  n steps ((xn-x0)/h-1); 
    double a;          // final y
 
    System.out.println("rk1_xy_acc.java running");
    System.out.println("yp(x,y) = 4*x^3 = y'");
    System.out.println("yf(x,y) = x^4");
    System.out.println("initial x0="+x0+"  y0="+y0+"  h="+h);
    a = rungekutta(x0, y0, xn, h); // 4th order
    System.out.println(" ");
    System.out.println(" ");
    a = third(x0, y0, xn, h); // third order
    System.out.println(" ");
    System.out.println(" ");
    a = second(x0, y0, xn, h); // second order
    System.out.println(" ");
    System.out.println(" ");
    a = first(x0, y0, xn, h); // first order
    System.out.println(" ");
    System.out.println(" ");
    System.out.println("rk1_xy_acc.java finished");
  } // end rk1_xy_acc

  double rungekutta(double x0, double y0, double xn, double h) // and yp(x,y)
  {
    int n = (int)((xn-x0)/h)-1;
    double k1, k2, k3, k4, k5;
    double y = y0;
    double x = x0;
    
    System.out.println("rungekutta running, n="+n+"  h="+h);
    System.out.println("x0="+x0+"  y0="+y0+"  chk="+chk(x0,y0));
    for(int i=0; i<=n; i++)
    {
      k1 = h*yp(x, y);
      k2 = h*yp(x+h/2.0, y+k1/2.0);
      k3 = h*yp(x+h/2.0, y+k2/2.0);
      k4 = h*yp(x+h, y+k3);
      k5 = (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
      System.out.println("k1="+k1+"  k2="+k2);
      System.out.println("k3="+k3+"  k4="+k4+"  k5="+k5);      
      y = y + k5;
      x = x + h;
      System.out.println("x="+x+"  y="+y+"  chk="+chk(x,y));
      System.out.println(" ");
    } // end i   
    return y;
  } // end rungekutta

  double first(double x0, double y0, double xn, double h) // and yp(x,y)
  {
    int n = (int)((xn-x0)/h)-1;
    double y = y0;
    double x = x0;
    
    System.out.println("first running, n="+n+"  h="+h);
    System.out.println("x0="+x0+"  y0="+y0+"  chk="+chk(x0,y0));
    for(int i=0; i<=n; i++)
    {
      y = y + h*yp(x,y);
      x = x + h;
      System.out.println("x="+x+"  y="+y+"  chk="+chk(x,y));
    } // end i   
    return y;
  } // end first

  double second(double x0, double y0, double xn, double h) // and yp(x,y)
  {
    int n = (int)((xn-x0)/h)-1;
    double k1, k2, k3;
    double y = y0;
    double x = x0;
    
    System.out.println("second running, n="+n+"  h="+h);
    System.out.println("x0="+x0+"  y0="+y0+"  chk="+chk(x0,y0));
    for(int i=0; i<=n; i++)
    {
      k1 = h*yp(x,y);
      k2 = h*yp(x+h,y+k1);
      k3 = (k1+k2)/2.0;
      y = y + k3;
      x = x + h;
      System.out.println("x="+x+"  y="+y+"  chk="+chk(x,y));
    } // end i   
    return y;
  } // end second

  double third(double x0, double y0, double xn, double h) // and yp(x,y)
  {
    int n = (int)((xn-x0)/h)-1;
    double k1, k2, k3, k4;
    double y = y0;
    double x = x0;
    
    System.out.println("third running, n="+n+"  h="+h);
    System.out.println("x0="+x0+"  y0="+y0+"  chk="+chk(x0,y0));
    for(int i=0; i<=n; i++)
    {
      k1 = h*yp(x,y);
      k2 = h*yp(x+h/2.0,y+k1);
      k3 = h*yp(x+h,y+k2);
      k4 = (k1+2.0*k2+k3)/4.0;
      y = y + k4;
      x = x + h;
      System.out.println("x="+x+"  y="+y+"  chk="+chk(x,y));
    } // end i   
    return y;
  } // end third

  double yp(double x, double y)
  { // return information on 1st derivative
      return 4.0*x*x*x;
  } // end yp

  double yf(double x, double y)
  { 
    return x*x*x*x;
  } // end yf

  double chk(double x, double y) // error
  { 
    return x*x*x*x - y;
  } // end chk

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