//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
// Basit diferansiyel denklemler dy/dx=f(x,y)

import java.io.*;

abstract class f_xy
{
  //single function single independent variable
  // example f=x*x
  abstract double func(double x,double y);
}

class fxy extends f_xy
{
  double func(double x,double y)
  {return -2*x*x*x+12*x*x-20*x+8.5;}
}

class fxy1 extends f_xy
{
  double func(double x,double y)
  {return 4.0*Math.exp(0.8*x)-0.5*y;}
}

class diferansiyelE
{
  public static double[][] euler(f_xy f,double x0,double y0,double xn,double h)
  {
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;
   a[1][i+1]=a[1][i]+f.func(a[0][i],a[1][i])*h;
   i++;
   }
  return a;  
  }
  
  public static double[][] heun(f_xy f,double x0,double y0,double xn,double h)
  {
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  double yi0=0;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;   
   yi0=a[1][i]+f.func(a[0][i],a[1][i])*h;
   a[1][i+1]=a[1][i]+(f.func(a[0][i],a[1][i])+f.func(a[0][i+1],yi0) )/2.0;
   i++;
   }
  return a;  
  }
  
  public static double[][] polygon(f_xy f,double x0,double y0,double xn,double h)
  {
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  double yi0=0;
  double xi0;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;
    xi0=xx+h/2.0;    
   yi0=a[1][i]+f.func(a[0][i],a[1][i])*h/2.0;
   a[1][i+1]=a[1][i]+f.func(xi0,yi0)*h;
   i++;
   }
  return a;  
  }
  
  public static double[][] RK2heun(f_xy f,double x0,double y0,double xn,double h)
  {
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  double k1,k2;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;
   k1=f.func(a[0][i],a[1][i]);
   k2=f.func(a[0][i]+h,a[1][i]+k1*h);   
   a[1][i+1]=a[1][i]+(k1+k2)/2.0*h;
   i++;
   }
  return a;  
  }
  
  public static double[][] RK2orta(f_xy f,double x0,double y0,double xn,double h)
  {
	  //orta nokta yöntemi
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  double k1,k2;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;
   k1=f.func(a[0][i],a[1][i]);
   k2=f.func(a[0][i]+0.5*h,a[1][i]+0.5*k1*h);   
   a[1][i+1]=a[1][i]+k2*h;
   i++;
   }
  return a;  
  }
  
  public static double[][] RK2(f_xy f,double x0,double y0,double xn,double h)
  {
	  //Ralston 2/3
  int n=(int)((xn-x0)/h);	
  System.out.println("n="+n);  
  double a[][]=new double[2][n+1];
  int i=0;
  a[0][0]=x0;
  a[1][0]=y0;
  double k1,k2;
  for(double xx=x0;xx<xn;xx+=h)
   {a[0][i+1]=xx+h;
   k1=f.func(a[0][i],a[1][i]);
   k2=f.func(a[0][i]+3.0/4.0*h,a[1][i]+3.0/4.0*k1*h);   
   a[1][i+1]=a[1][i]+(1.0/3.0*k1+2.0/3.0*k2)*h;
   i++;
   }
  return a;  
  }
  
  
  public static void main(String args[]) throws IOException
  {double x0=0;
   double y0=1.0;
   fxy f1=new fxy();
   fxyKC fKC=new fxyKC();
   double a[][]=RK2(f1,x0,y0,4.0,1.0);
   Text.print(Text.T(a));
   Plot p1=new Plot(a);
   p1.setPlotType(0,24);
   p1.addSpline(a,8);
   p1.addFunction(fKC,0.0,4.0,200);
   p1.plot();
  }
}