// Continuation-homotopy
// roots of nonlinear system of equations
// Referance : Numerical Analysis, 8th edition Richard L. Burden, J. Douglas Faires 3
// page 637
// All the required derivatives calculated numercally
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Function definition ================

class fa extends fi_xi
{
public double[] func(double x[])
{	
//çözümü istenen fonksiyon	
double ff[]=new double[3];
ff[0]=3.0*x[0]-Math.cos(x[1]*x[2])-0.5;
ff[1]=x[0]*x[0]-81.0*(x[1]+0.1)*(x[1]+0.1)+Math.sin(x[2])+1.06;
ff[2]=Math.exp(-x[0]*x[1])+20.0*x[2]+(10.0*Math.PI-3.0)/3.0;
return ff;
} }


class fb extends fij_xi
{
fi_xi f;	

public fb(fi_xi fi)
{f=fi;}

public double[][] func(double x[])
{ 
  int n=x.length;
  double a[][]=new double[n][n];
  for(int i=0;i<n;i++)
  {
    for(int j=0;j<n;j++)
    {
        a[i][j]=derivative(x,i,j);
    }
  }
return a;
} 

public double derivative(double x[],int denklem_ref,int x_ref)
{
// fonksiyon f in denklem_ref sayılı fonksiyonunun(sayılar 0 dan başlar)
// x[x_ref} değişkenine göre türevi (sayılar 0 dan başlar) 
// df_denklem_ref(x)/d_x_ref
// bu metod newtond metodu içinde kullanılmak içindir.
double h0=0.0256;
int i,m;
int n=7;
double f1[];
f1=new double[x.length];
double f2[];
f2=new double[x.length];
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
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++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
f1=f.func(x1);
f2=f.func(x2);
T[i][0]=( f1[denklem_ref] - f2[denklem_ref])/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
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 class NA30B
{
// Süreklilik(Continuity-homotopy) metodu kullanarak denklem sisteminin köklerini bulmak	
// Bu metod denklem sistemi çözümü için Matrix.java programını
// kullanmaktadır.


public static double[] multiply(double left,double[] right)
{
//multiplying a vector with a constant
int i;
int n=right.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
  {b[i]=left*right[i];}
return b;
}

public static double[] multiply(double[][] left,double[] right)
{
//multiplication of one matrix with one vector
int ii,jj,i,j,k;
int m1=left[0].length;
int n1=left.length;
int m2=right.length;
double[] b;
b=new double[m2];
 if(n1 != m2)
 {
 System.out.println("inner matrix dimensions must agree");
 for(ii=0;ii<n1;ii++)
   {
     b[ii]=0;
   }
   return b;
 }
 for(i=0;i<m1;i++)
   {
     b[i]=0;
     for(k=0;k<n1;k++)
	b[i]+=left[i][k]*right[k];
   }
return b;
//end of multiply of a matrix and a vector
}

public static double[] add(double[] left,double[] right)
{
//addition of two vectors
int n1=left.length;
int n2=right.length;
int nMax;
int i;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[];
b=new double[nMax];
for(i=0;i<n1;i++)
  {
    b[i]=b[i]+left[i];
  }
for(i=0;i<n2;i++)
  {
    b[i]=b[i]+right[i];
  }
return b;
//end of vector addition method
}

public static double[][] inv(double[][] a)
{
// INVERSION OF A MATRIX
// inversion by using gaussian elimination
// with full pivoting
int n=a.length;
int m=a[0].length;
double b[][];
b=new double[n][n];
int indxc[];
int indxr[];
double ipiv[];
indxc=new int[n];
indxr=new int[n];
ipiv=new double[n];
int i,j,k,l,ll,ii,jj;
int icol=0;
int irow=0;
double big,dum,pivinv,temp;
if(n!=m)
{
   System.out.println("Matrix must be square ");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
	  return b;
}
for(i=0;i<n;i++)
   for(j=0;j<n;j++)
       b[i][j]=a[i][j];
for(i=0;i<n;i++)
{
  big=0.0;
  for(j=0;j<n;j++)
  {
    if(ipiv[j] != 1)
      for(k=0;k<n;k++)
      {
	if(ipiv[k] == 0)
	{
	  if(Math.abs(b[j][k]) >= big)
	  {
	    big=Math.abs(b[j][k]);
	    irow=j;
	    icol=k;
	  }
	}
	else if(ipiv[k] > 1 )
	  {
	  System.out.println("error : inverse of the matrix : singular matrix-1");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
	  return b;
	  }
      }
  }
  ++ ipiv[icol];
  if(irow != icol)
    for(l=0;l<n;l++)
    {
    temp=b[irow][l];
    b[irow][l]=b[icol][l];
    b[icol][l]=temp;
    }
  indxr[i]=irow;
  indxc[i]=icol;
  if(b[icol][icol] == 0.0)
  {
    System.out.println("error : inverse of the matrix : singular matrix-2");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
    return b;
  }
  pivinv=1.0/b[icol][icol];
  b[icol][icol]=1.0;
  for(l=0;l<n;l++) b[icol][l] *=pivinv;
  for(ll=0;ll<n;ll++)
  if(ll != icol)
  {
    dum=b[ll][icol];
    b[ll][icol]=0.0;
    for(l=0;l<n;l++) b[ll][l]-= b[icol][l]*dum;
  }
}
for(l=n-1;l>=0;l--)
{
  if(indxr[l] != indxc[l])
     for(k=0;k<n;k++)
     {
       temp=b[k][indxc[l]];
       b[k][indxc[l]]=b[k][indxr[l]];
       b[k][indxr[l]]=temp;
     }
}
return b;
}
  
public static double[] continuationRK4(fi_xi f,fij_xi J,double x[],int N)
{
//====================================================	
// Roots of nonlinear system of equations Homotophy RK6
// 4.th order Runge-Kutta 
// yi+1 = yi + (1/6)*( k1 + 2k2 +2k3+k4)h
// k1=f(xi,yi)
// k2=f(xi+0.5h,yi+0.5k1h)
// k3=f(xi+0.5h,yi+0.5k2h)
// k4=f(xi+h,yi +k3h)
//===================================================
//x vector of independent variables
//y vector of dependent variables
//dy derivative vector of dependent variables
int i;
int nmax=400;
double tolerance=1.0e-20;
int n=x.length;
double h=1.0/(double)N;
double b[]=new double[n];
double x1[]=new double[n];
double k[][]=new double[4][n];
double A[][]=new double[n][n];
b=multiply(-h,f.func(x));
for(i=0;i<N;i++)
{
x1=x;	
A=J.func(x1);
k[0]=multiply(inv(A),b);
x1=add(x,multiply(0.5,k[0]));
A=J.func(x1);
k[1]=multiply(inv(A),b);
x1=add(x,multiply(0.5,k[1]));
A=J.func(x1);
k[2]=multiply(inv(A),b);
x1=add(x,k[2]);
A=J.func(x1);
k[3]=multiply(inv(A),b);
for(int j=0;j<n;j++) {x[j]=x[j]+1.0/6.0*(k[0][j]+2.0*k[1][j]+2.0*k[2][j]+k[3][j]);}
}
return x;
}

public static double[] continuationRK6(fi_xi f,fij_xi J,double x[],int N)
{
//====================================================	
// Roots of nonlinear system of equations Homotophy RK6
// yi+1 = yi + (1/90)*( 7k1 + 32k3 +12k4+32k5+7k6)h
// k1=f(xi,yi)
// k2=f(xi+0.25h , yi+0.25k1h)
// k3=f(xi+0.25h , yi+0.125k1h+0.125k2h)
// k4=f(xi+0.5h , yi - 0.5k2h+k3h)
// k5=f(xi+0.75h , yi + (3/16)k1h+(9/16)k4h)
// k6=f(xi+h , yi - (3/7)k1h+(2/7)k2h+(12/7)k3h - (12/7)k4h+(8/7)k5h)
//===================================================
//x vector of independent variables
//y vector of dependent variables
//dy derivative vector of dependent variables
int i;
int nmax=400;
double tolerance=1.0e-20;
int n=x.length;
double h=1.0/(double)N;
double b[]=new double[n];
double x1[]=new double[n];
double k[][]=new double[6][n];
double A[][]=new double[n][n];
b=multiply(-h,f.func(x));
for(i=0;i<N;i++)
{
A=J.func(x);
// k1=f(xi,yi)
k[0]=multiply(inv(A),b);
x1=add(x,multiply(0.25,k[0]));
A=J.func(x1);
// k2=f(xi+0.25h , yi+0.25k1h)
k[1]=multiply(inv(A),b);
x1=add(x,add(multiply(0.125,k[0]),multiply(0.125,k[1])));
A=J.func(x1);
// k3=f(xi+0.25h , yi+0.125k1h+0.125k2h)
k[2]=multiply(inv(A),b);
x1=add(x,add(multiply(-0.5,k[1]),k[2]));
A=J.func(x1);
// k4=f(xi+0.5h , yi - 0.5k2h+k3h)
k[3]=multiply(inv(A),b);
x1=add(x,add(multiply((3.0/16.0),k[0]),multiply((9.0/16.0),k[3])));
A=J.func(x1);
// k5=f(xi+0.75h , yi + (3/16)k1h+(9/16)k4h)
k[4]=multiply(inv(A),b);
x1=add(x,
       add(multiply((-3.0/7.0),k[0]),add(multiply((2.0/7.0),k[1]),
       add(multiply((12.0/7.0),k[2]),
           add(multiply((-12.0/7.0),k[3]),multiply((8.0/7.0),k[4]))))));
A=J.func(x1);
// k6=f(xi+h , yi - (3/7)k1h+(2/7)k2h+(12/7)k3h - (12/7)k4h+(8/7)k5h)
k[5]=multiply(inv(A),b);
// yi+1 = yi + (1/90)*( 7k1 + 32k3 +12k4+32k5+7k6)h
for(int j=0;j<n;j++) {x[j]=x[j]+1.0/90.0*(7.0*k[0][j]+32.0*k[2][j]+12.0*k[3][j]+32.0*k[4][j]+7.0*k[5][j]);}
}
return x;
}
public static double[] newton_continuationRK4(fi_xi f,fij_xi df,double x[])
{
// Roots of nonlinear system of equations Homotophy RK4 plus Newton-Raphson
// ti : weight function
// x vector of independent variables
// y vector of dependent variables
// dy derivative vector of dependent variables
x=continuationRK4(f,df,x,4);
double ti=1.0;
int i;
int nmax=400;
double tolerance=1.0e-10;
int n=x.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
{
b[i]=1.0;
}
i=0;
while( i++ < nmax && Matrix.abs(b) > tolerance )
{ b=Matrix.multiply(Matrix.divide(f.func(x),df.func(x)),-ti); 
  x=Matrix.add(x,b);
}
if(i >= nmax) JOptionPane.showMessageDialog(null,"Warning : Maximum number of iterations are exceeded \n"+
              " Results may not be valid","MAXIMUM ITERATIONS WARNING",JOptionPane.WARNING_MESSAGE);
System.out.println();
return x;
}
public static double[] enterdata()
{
    String s=JOptionPane.showInputDialog("enter initial root estimations with a space in between ");        
    StringTokenizer token=new StringTokenizer(s);
    int n=token.countTokens()-1;
    int m=n+1;
    double a[]=new double[m];
    int j=0;           
    while(token.hasMoreTokens())
    {
    Double ax=new Double(token.nextToken());
    a[j++]=ax.doubleValue();
    }
    return a; 
}
     public static void main (String args[])  
     {
	 fa f =new fa();
     fb df=new fb(f);
     double [] x0=enterdata();
     double [] r= continuationRK4(f,df,x0,4);
     double [] r1= continuationRK6(f,df,x0,4);
     double [] r2= newton_continuationRK4(f,df,x0);
     String s=" Root values with homotophy RK4 : \n"+Matrix.toStringT(r);
     s+=" Root values with homotophy RK6 : \n"+Matrix.toStringT(r1);
     s+=" Root values with homotophy RK4 plus Newton-Raphson : \n"+Matrix.toStringT(r2);
     String s1="Finding multiple roots with homotophy method : ";
     JOptionPane.showMessageDialog(null,s,s1,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}
