// SCO5B Süreklilik(Continuation-homotopy) metoduyla
// birden fazla değişkenli denklem sisteminin köklerini bulmak
// Metod Referansı : Numerical Analysis, 8th edition Richard L. Burden, J. Douglas Faires 3
// sayfa 637
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Tanımı gereken fonksiyon ================
abstract class fi_xi
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x) returns the value of f[0] and f[1]
  // as a two dimensional vector

  abstract double[] func(double x[]);
}

abstract class fij_xi
{
  // multifunction multi independent variable
  // for n independent variable n*n matrix of
  // dependent variables are returned
  // example f[0][0]=x[0]+Math.sin(x[1])   f[0][1]=x[0]-x[1]
  //         f[1][0]=x[0]*x[0]-x[1]        f[1][1]=Math.exp(x[0]+x[1]*x[1]
  // func(x) returns the value of f[0][0], f[0][1]
  //                              f[1][0], f[1][1]
  // as a two dimensional matrix
  abstract double[][] func(double x[]);
}

class fa extends fi_xi
{
public double[] func(double x[])
{
//çözümü istenen fonksiyonun türev fonksiyonları	
//çö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
{
public double[][] func(double x[])
{
//çözümü istenen fonksiyonun türev fonksiyonları	
double a[][]=new double[3][3];
a[0][0]=3.0;     a[0][1]=x[2]*Math.sin(x[1]*x[2]);a[0][2]=x[1]*Math.sin(x[1]*x[2]);
a[1][0]=2.0*x[0];a[1][1]=-162.0*(x[1]+0.1);a[1][2]=Math.cos(x[2]);
a[2][0]=-x[1]*Math.exp(-x[0]*x[1]);a[2][1]=-x[0]*Math.exp(-x[0]*x[1]);a[2][2]=20.0;
return a;
} 
}

public class FSCO5B
{
// 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)
{
//====================================================	
// lineer olmıyan denklem siteminin kökleri
// fj(xi)=0 lineer olmıyan denklem sistemi için çözüm matrisi 
// türev matrisi J de tanımlanacaktır
// Homotopy, continuation (devamlılık) metodu ile lineer olmıyan denklem sistemi çözümü 
// 4.derece Runge-Kutta denklemi ile
// 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  : bağımsız değişken vektörü
// y  : bağımlı değişken vektörü
// dy : bağımlı değişkenin türev fonksiyonları vektörü
// N  : Runge-Kutta diferansiyel denklem noktasayısı

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++)
{
A=J.func(x);
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)
{
//====================================================	
// lineer olmıyan denklem siteminin kökleri
// fj(xi)=0 lineer olmıyan denklem sistemi için çözüm matrisi 
// türev matrisi J de tanımlanacaktır
// Homotopy, continuation (devamlılık) metodu ile lineer olmıyan denklem sistemi çözümü 
// 6. derece Runge-Kutta denklemi ile
// 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  : bağımsız değişken vektörü
// y  : bağımlı değişken vektörü
// dy : bağımlı değişkenin türev fonksiyonları vektörü
// N  : Runge-Kutta diferansiyel denklem noktasayısı
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[])
{
// lineer olmıyan denklem siteminin kökleri
// fj(xi)=0 lineer olmıyan denklem sistemi için çözüm matrisi 
// türev matrisi dfj/dxi de tanımlanacaktır
// Newton-Raphson metodu ile lineer olmıyan denklem sistemi çözümü 
//ti :önem katsayısı
//x bağımsız değişken vektörü
//y bağımlı değişken vektörü
//dy bağımlı değişkenin türev fonksiyonları vektörü
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,"Uyarı maximum iterasyon sayısı aşıldı \n"+
              " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE);
System.out.println();
return x;
}

public static double[] verigir()
{
    String s=JOptionPane.showInputDialog("kök tahmin değerlerini arada bir boşluk bırakarak giriniz: ");        
    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();
     double [] x0=verigir();
     double [] r= continuationRK4(f,df,x0,4);
     double [] r1= continuationRK6(f,df,x0,4);
     double [] r2= newton_continuationRK4(f,df,x0);
     String s=" kök değerleri süreklilik RK4 : \n"+Matrix.toStringT(r);
     s+=" kök değerleri süreklilik RK6 : \n"+Matrix.toStringT(r1);
     s+=" kök değerleri newton-süreklilik RK4 : \n"+Matrix.toStringT(r2);
     String s1="çok değişkenli süreklilik metodu ile kök bulma : ";
     JOptionPane.showMessageDialog(null,s,s1,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}