//======================================================
// Numerical Analysis package in java
// example to show utilisation of integration (integral)
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;


class f1 extends f_x
{ public double func(double x)
{return x*x*x;}
}
class dif1
{
	
public static double derivative(f_x f,double x)
{
double h0=0.001;
int i,m;
int n=7;
double T[][]=new double[n][n];
double h[];
h=new double[n];
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++) {h[i]=h0*Math.pow(r,i);}
for(i=0;i<n;i++)  {T[i][0]=dfdx(f,x,h[i]);}
for(m=1;m<n;m++)
{ for(i=0;i<n-m;i++)
  {T[i][m]=(h[i]*h[i]*T[i+1][m-1]-h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]- h[i+m]*h[i+m]);}
}
double xx=T[0][n-1];
return xx;
}

public static double dfdx_ridder(f_x f, double x, double h)
{
double CON=1.4;
double CON2=CON*CON;
double BIG=1.0e30;
int NTAB=10; 
double SAFE=2.0;
int i,j;
double errt,fac,hh;
double ans=0;
double hata;
if (h == 0.0) System.out.println("h must be nonzero in dfridr.");
double a[][]=new double[NTAB][NTAB];
hh=h;
a[1][1]=(f.func(x+hh)-f.func(x-hh))/(2.0*hh);
hata=BIG;
for (i=2;i<=NTAB;i++) 
{
hh /= CON;
a[1][i]=(f.func(x+hh)-f.func(x-hh))/(2.0*hh);
fac=CON2;
for (j=2;j<=i;j++) { 
a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
fac=CON2*fac;
errt=Math.max(Math.abs(a[j][i]-a[j-1][i]),Math.abs(a[j][i]-a[j-1][i-1]));
if (errt <= hata) 
{ 
hata=errt;
ans=a[j][i];
}
}
if (Math.abs(a[i][i]-a[i-1][i-1]) >= SAFE*hata) break;
}
return ans;
}

public static double dfdx_richardson(f_x f,double x)
{
double h0=0.01;
double h=h0;
int i,m;
int n=10;
double T[][]=new double[n][n];
for(i=0;i<n;i++)  {T[i][0]=dfdx(f,x,h);h/=2.0;}
for(m=1;m<n;m++)
{ for(i=0;i<n-m;i++)
  {T[i][m]=(4.0/3.0*T[i+1][m-1]-1.0/3.0*T[i][m-1]);}
}
double xx=T[0][n-1];
return xx;
}


public static double dfdx(f_x f,double x)
{ double h=0.00001;
  return dfdx(f,x,h);
}
	
public static double dfdx(f_x f,double x,double h)
{ double hh=1.0/h; 
  return (-f.func(x+2.0*h)+8.0*f.func(x+h)-8.0*f.func(x-h)+f.func(x-2.0*h))/12.0*hh;
}

public static double d2fdx2(f_x f,double x)
{ double h=0.00001;
  return d2fdx2(f,x,h);
}

public static double d2fdx2(f_x f,double x,double h)
{
  double hh=1.0/h; 
	return (-f.func(x+2.0*h)+16.0*f.func(x+h)-30.0*f.func(x)+16.0*f.func(x-h)-f.func(x-2.0*h))/(12.0)*hh*hh;}

public static double d3fdx3(f_x f,double x)
{ double h=0.00001;
  return d3fdx3(f,x,h);}

public static double d3fdx3(f_x f,double x,double h)
{  double hh=1.0/h;
	return (-f.func(x+3.0*h)+8.0*f.func(x+2.0*h)-13.0*f.func(x+h)+13.0*f.func(x-h)-8.0*f.func(x-2.0*h)+f.func(x-3.0*h))/(8.0)*hh*hh*hh;}

  public static void main(String args[]) 
  {
  //Gauss-Legendre integrali
  f1 b1=new f1();
  System.out.println("turev dfdx    : "+dfdx(b1,1.0));
  System.out.println("turev dfdx (ridder formülü)    : "+dfdx_ridder(b1,1.0,0.005));
  System.out.println("turev d2dfx2  : "+d2fdx2(b1,1.0));
  System.out.println("turev d3dfx3  : "+d3fdx3(b1,1.0));
  
  System.out.println("turev dfdx (romberg formülü)   : "+derivative(b1,1.0));
  System.out.println("turev dfdx (richardson formülü)   : "+dfdx_richardson(b1,1.0));
  }
}
