//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
//Blasius flat plate problem

import java.io.*;

class fm1 extends f_xi
{
//multivariable function
 double func(double x[],int x_ref)
 {
   //Blasius flat plate diferential equation
   //f'''+0.5*ff'' = 0
   //f'''=-0.5*x[1]*x[3] f''=x[3] f'=x[2] f=x[1]]
   double a=0;
   if(x_ref==2) a= -0.5*x[1]*x[3];
   if(x_ref==1) a= x[3];
   if(x_ref==0) a= x[2];
   return a;
 }
}


class diferansiyel1A
{

public static double[][] RK4(f_xi fp,double x0,double xn,double f0[],int N)
{
//4th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[4];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[j]=fi[j-1][i]+k[1]/2.0;
     k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]+k[2];
     k[3]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][]
RKF45(f_xi fp,double x0,double xn,double f0,int N) throws IOException
{
// Runge Kutta - Fehlberg Method
// fp  : derivative function df/dx(f,x)  (defined as a class)
// Tol :Tolerance
//Hmax : maximum possible step size
//Hmin : minimum possible step size
//k    : Runge kutta coefficients
//Err  : error
//x[i] : input variable to fp = df/dx(f,x) = df/dx(x[i])
//j    : actual step size
//fi[][]:solution matrix
//fj[][]:solution matrix in exact size(j) (return function)
double h=(xn-x0)/N;
double fi[][]=new double[2][1000];
double Tol=2.0e-3;
double Hmin=h/64.0;
double Hmax=h*64.0;
double k[]=new double[6];
double Err;
double s;
double xi[]=new double[2];
int j=0;
fi[0][j]=x0;
fi[1][j]=f0;
while(fi[0][j] < xn )
  {
  if(fi[0][j]+h > xn) h=xn-fi[0][j];
  xi[0]=fi[0][j];
  xi[1]=fi[1][j];
   k[0]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h/4.0;
  xi[1]=fi[1][j]+k[0]/4.0;
   k[1]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+3.0/8.0*h;
  xi[1]=fi[1][j]+3.0/32.0*k[0]+9.0/32.0*k[1];
   k[2]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+12.0/13.0*h;
  xi[1]=fi[1][j]+1932.0/2197.0*k[0]-7200.0/2197.0*k[1]+7296.0/2197.0*k[2];
   k[3]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h;
  xi[1]=fi[1][j]+439.0/216.0*k[0]-8.0*k[1]+3680.0/513.0*k[2]-845.0/4104.0*k[3];
   k[4]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+0.5*h;
  xi[1]=fi[1][j]-8.0/27.0*k[0]+2.0*k[1]
  -3544/2565*k[2]+1859.0/4104.0*k[3]-11.0/40.0*k[4];
   k[5]=h*fp.func(xi,0);
   Err=Math.abs(1.0/360.0*k[0]-128/4275*k[2]
   -2197.0/75240*k[3]+1.0/5.0*k[4]+2.0/55.0*k[5]);
  if(Err<Tol || h<2*Hmin)
     {
     fi[0][j+1]=fi[0][j]+h;
     fi[1][j+1]=fi[1][j]+16.0/135.0*k[0]+6656.0/12825.0*k[2]+
     28561.0/56430.0*k[3]-9.0/50.0*k[4]+2.0/55.0*k[5];
     j++;
     }
  else
     {
     if(Err==0) s=0;
     else s=0.84*(Math.pow(Tol*h/Err,0.25));
     if(s<0.75 && h>(2*Hmin)  ) h/=2.0;
     if(s>1.5  && Hmax >(2*h) ) h*=2.0;
     }
  }
double fj[][]=new double[2][j+1];
for(int i=0;i<=j;i++)
  {
  fj[0][i]=fi[0][i];
  fj[1][i]=fi[1][i];
  }
return fj;
}

  public static double [] thomas(double a[][],double r[])
  {
  //
  int n=a.length;
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double x[]=new double[n];
  for(int i=0;i<n;i++)     {f[i]=a[i][i];}
  for(int i=0;i<(n-1);i++) {g[i]=a[i][i+1];}
  for(int i=0;i<(n-1);i++) {e[i+1]=a[i+1][i];}
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  public static double [] thomas(double f[],double e[],double g[],double r[])
  {
  int n=f.length;
  double x[]=new double[n];
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  
  public static double [][] cubic_spline(double xi[],double yi[],double c0,double cn)
  {
  int n=xi.length;
  double h[]=new double[n];
  double w[]=new double[n];
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double d[]=new double[n];
  double x[]=new double[n];
  double S[][]=new double[4][n];
  int k;	  
  for(k=0;k<(n-1);k++)
     {h[k]=xi[k+1]-xi[k];
      w[k]=(yi[k+1]-yi[k])/h[k];      
	 }
  d[0]=c0;
  d[n-1]=cn;
  for(k=1;k<(n-1);k++)
     {d[k]=6.0*(w[k]-w[k-1]);}
  f[0]=1.0;
  f[n-1]=1.0;
  g[0]=0.0;
  g[n-1]=0.0;
  e[0]=0.0;
  e[n-1]=0.0;
  for(k=1;k<(n-1);k++)
  {f[k]=2.0*(h[k]+h[k-1]);e[k]=h[k-1];g[k]=h[k];}
  S[2]=thomas(f,e,g,d);
  S[3]=xi;
  for(k=0;k<(n-1);k++)
    {S[0][k]=(6.*yi[k+1]-h[k]*h[k]*S[2][k+1])/(6.0*h[k]);
     S[1][k]=(6.*yi[k]-h[k]*h[k]*S[2][k])/(6.0*h[k]);
    } 
  return S;    	 	     	 
  }
  
  public static double funcSpline(double S[][],double x)
  {
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  }     
   
    public static double[][]  funcSpline(double xi[],double yi[],int aradegersayisi)
    {
	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double S[][]=cubic_spline(xi,yi,0,0);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {  z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);k++;
          for(int j=0;j<aradegersayisi;j++) 
             {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
          z[0][k]=z[0][k-1]+dx;z[1][k]=funcSpline(S,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);
    return z;
    } 
 	
  public static void main(String args[]) throws IOException
  {
  fm1 b3=new fm1();
  double y[]=new double[3];
  y[0]=0.0; 
  y[1]=0.0; 
  y[2]=0.3; //y0 first guess for f''
  // solution of the differential system : 
  double xi[]=new double[30];
  double fi[]=new double[30];
  double a[][];
  // shooting method
  for(int i=1;i<30;i++)
  {
  a=RK6(b3,0.0,6.8,y,10000);
  xi[i]=a[2][10000];
  fi[i]=y[2];
  y[2]+=0.1;  //change the boundary value
  }
  //curve fit to find actual boundary value
  double S[][]=cubic_spline(xi,fi,0,0);
  y[2]=funcSpline(S,1.0);
  //y[0]=Numeric.f_least_square(Numeric.poly_least_square(xi,fi,4),1.0);//curve fit value for f''
  // calculate final solution
  a=RK6(b3,0.0,6.8,y,10000);
  String heading[]={"eta","fi","fi'","fi''"};
  Text.print(Text.T(a),heading);
  Plot pp=new Plot(a);
  pp.plot();
  }
}
