// FSCO5H Newton Raphson metoduyla
// birden fazla değişkenli denklem sisteminin köklerini bulmak
// ve optimizasyon Hessian matris H sayısal olarak hesaplanmış

import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Tanımı gereken fonksiyon ================
abstract class fij_xi extends Mathd
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0][0]=x[0]+sin(x[1])  f[0][1]=x[0]*x[1]
  //         f[1][0]=x[0]*x[0]-x[1]  f[1][1]=x[0]+x[1]*x[0]
  // func(x) returns the value of f[0][0] and f[0][1]
  //                              f[0][1] and f[1][1]
  // as a two dimensional vector
  //
  abstract double[][] func(double x[]);
}

abstract class fi_xi extends Mathd
{
  // 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][0] 
  //                              f[0][1] 
  // as a two dimensional vector
  //
  abstract double[] func(double x[]);
}

abstract class f_xj extends Mathd
{
  // single function  multi independent variable
  // a single value is returned indiced to equation_ref
  // example f=x[0]+sin(x[1])
  //         
  // func(x,1) returns the value of f[1]
  // func(x,0) returns the value of f[0]

  abstract double func(double x[]);
}

class fxa extends f_xj
{
public double func(double x[]) 
{
//çözümü istenen fonksiyon	
double ff;
ff=3.0*x[0]*x[0]-4.0*x[0]*x[1]+2.0*x[1]*x[1]-x[0]-x[1];
return ff; //minimum değil maxsimum istendiğinden - ile çarptık
}
}

class fxb extends fi_xi
{
f_xj fxx;	

fxb(f_xj fxi) 
{fxx=fxi;}
 	
public double[] func(double x[])
{
// Türev hesabı
double h=0.0013421432165;
double hh=1.0/h;
int i,m;
int n=x.length;
double f1;
double f2;
double f3;
double f4;
double x1[];
x1=new double[n];
double x2[];
x2=new double[n];
double x3[];
x3=new double[n];
double x4[];
x4=new double[n];
double fx[]=new double[n];

for(int j=0;j<n;j++)
{
for(i=0;i<n;i++)
{
x1[i]=x[i];
x2[i]=x[i];
x3[i]=x[i];
x4[i]=x[i];
}
x1[j]+=2.0*h;
x2[j]+=h;
x3[j]-=h;
x4[j]-=2.0*h;

f1=fxx.func(x1);
f2=fxx.func(x2);
f3=fxx.func(x3);
f4=fxx.func(x4);
fx[j]=(-f1+8.0*f2-8.0*f3+f4)/12.0*hh;
}
return fx;
} 

}

public class FSCO5J
{
// Newton Raphson Metodunu kullanarak denklem sisteminin köklerini bulmak	
// Bu metod denklem sistemi çözümü için Matrix.java programını
// kullanmaktadır.


public static double[] newtond( f_xj f,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 turev alt programıyla hesaplanı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ü
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-15;
int n=x.length;
double b[];
b=new double[n];
double dy[][];
dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    b=multiply(divide(dy(f,x),H(f,x)),-ti);
    x=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);
return x;
}

public static double[] dy(f_xj fxx,double x[])
{
// Türev hesabı
double h=0.0013421432165;
double hh=1.0/h;
int i,m;
int n=x.length;
double f1;
double f2;
double f3;
double f4;
double x1[];
x1=new double[n];
double x2[];
x2=new double[n];
double x3[];
x3=new double[n];
double x4[];
x4=new double[n];
double fx[]=new double[n];

for(int j=0;j<n;j++)
{
for(i=0;i<n;i++)
{
x1[i]=x[i];
x2[i]=x[i];
x3[i]=x[i];
x4[i]=x[i];
}
x1[j]+=2.0*h;
x2[j]+=h;
x3[j]-=h;
x4[j]-=2.0*h;

f1=fxx.func(x1);
f2=fxx.func(x2);
f3=fxx.func(x3);
f4=fxx.func(x4);
fx[j]=(-f1+8.0*f2-8.0*f3+f4)/12.0*hh;
}
System.out.println(Matrix.toStringT(fx));
return fx;
}
 
public static double dfx(f_xj f,double x[],int denklem_ref,int x_ref)
{ // Türev hesabı
double h=0.000014354523465;
double hh=1.0/h;
int i,m;
int n=x.length;
double f1;
double f2;
double f3;
double f4;
double f5;

double x1[];
x1=new double[n];
double x2[];
x2=new double[n];
double x3[];
x3=new double[n];
double x4[];
x4=new double[n];
double x5[];
x5=new double[n];

for(i=0;i<n;i++)
{
x1[i]=x[i];
x2[i]=x[i];
x3[i]=x[i];
x4[i]=x[i];
x5[i]=x[i];
}
x1[x_ref]+=2.0*h;
x2[x_ref]+=h;
x4[x_ref]-=h;
x5[x_ref]-=2.0*h;

f1=f.func(x1);
f2=f.func(x2);
f3=f.func(x3);
f4=f.func(x4);
f5=f.func(x5);
double fx=(-f1+16.0*f2-30.0*f3+16.0*f3-f4)/12.0*hh*hh;
//System.out.println("f("+x_ref+","+denklem_ref+") = "+fx);
return fx;
}

public static double[][] H(f_xj f,double x[])
{
  int n=x.length;	
  double dy[][]=new double[n][n];;
  for(int ii=0;ii<n;ii++)
  {
    for(int jj=0;jj<n;jj++)
    {   //Hessian matrisi
        dy[ii][jj]=dfx(f,x,ii,jj);
    }
  }
  System.out.println(Matrix.toString(dy));
  return dy;
}

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


public static double[] divide(double[] left,double[][] right)
{
return pivotlugauss(right,left);
}
  //en kucuk kareler metodu	
  public static double[] pivotlugauss(double a[][],double b[])
  { //kısmi pivotlu gauss eleme yöntemi
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss eleme
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivotlama
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	}
	//gauss elemeyi çözme  
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  //geriye doğru yerine koyma
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
}

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[] 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[])  
     {
	 fxa f1 =new fxa();
     double [] x0=verigir();
     String s2="Newton-Raphson kök bulma : ";
     double [] r1= newtond(f1,x0);
     String s=" kök değerleri  : \n"+Matrix.toStringT(r1);
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}