BÖLÜM 12: ÖRNEKLERLE SAYISAL ANALİZ
12.1 SAYISAL ANALİZE GİRİŞ
Bu bölümümüzde bazı pratik matematik problemlerinin java ile nasıl çözülebileceğine göz atacağız. Sayısal analiz işlemlerini yapmak amacıyla iki sınıf oluşturduk bunlar Matrix ve Numeric sınıflarıdır. Bu serinin devamı olan java ile numerik analiz kitabında bu sınıflardaki metotları tek tek incelemeye çalışacağım. Burada ise sadece program listeleri ve örnek problemler vermekle yetineceğiz. Buradaki metodlar uluslararası projelerde kullanıldığından metot isim ve tanımlamaları ingilizcedir. Sadece bu kitabın kapsamı için değil gerçek nümerik problemlerde kullanılmak amacıyla geliştirilmişlerdir. Burada ayrıca hemen şunu da belirtelim: sayısal metodların kullanımı bu metodlar, sayısal analiz ve bilgisayar değişkenlerinin değişimi ile ilgili iyi bilgi gerektirmektedir. Bu sınıf ve metodları ciddi bir uygulamada kullanmadan önce verilen metodlar hakkında nümerik analiz kitaplarndan yeterince bilgi sahibi olmanızı ve test fonksiyonları kullanarak sizin istediğiniz görevi yapıp yapamıyacaklarını iyice kontrol etmenizi tavsiye ederim. Türkiyedeki nümerik analizle ilgilenen arkadaşlar eğer ilgilenirlerse benim kendi projelerimde kullanmak için geliştirdiğim bu gurubu daha da geliştirerek genel bir nümerik analiz kütüphanesine dönüştürmemiz mümkün olabilir.
12.2 MATRIX SINIFI
Matrix sınıfı double ve complex vektör (bir boyutlu değişken) ve matris (iki boyutlu değişken) değişkenlerinin matematik işlemlerini yapmak amacıyla tanımlanmıştır. complex (kompleks değişken sınıfı) daha önce türkçe olarak tanımlanan kompleks sınıfının aynıdır, fakat metot isimleri burada ingilizce olduğundan bu sınıfın tanımını tekrar verelim :
Program 12.1 complex.java, ingilizce metot isimleriyle kompleks değişkenler sınıfı .
//===========================================
// Numerical Analysis Package in Java
// Complex class defination
// Dr. Turhan Coban
//===========================================
import java.io.*;
// class complex
// complex number abstraction
//
class complex {
// constructors
public complex()
{
nreal=0;
nimag=0;
}
public complex(double nre,double nim)
{
nreal=nre;
nimag=nim;
}
public complex(double numer)
{
nreal=numer;
nimag=0;
}
public complex(complex value )
{
nreal=value.real();
nimag=value.imaginary();
}
// accessor functions
public double real()
{
return nreal;
}
public double imaginary()
{
return nimag;
}
public void setReal(double r)
{
nreal=r;
}
public void setImaginary(double i)
{
nimag=i;
}
public double R()
{
return Math.sqrt(nreal*nreal+nimag*nimag);
}
public double theta()
{
return Math.atan2(nimag,nreal);
}
public double dtheta()
{
return Math.atan2(nimag,nreal)*45.0/Math.atan(1.0);
}
// assignments
public void assign(complex right)
{
nreal=right.real();
nimag=right.imaginary();
}
public void assign(double nr,double ni)
{
nreal=nr;
nimag=ni;
}
public void add(complex right)
{
nimag = nimag + right.imaginary();
nreal = nreal + right.real();
}
public void substract(complex right)
{
nimag = nimag - right.imaginary();
nreal = nreal - right.real();
}
public void multiply(complex right )
{
nreal = nreal*right.real() - nimag*right.imaginary();
nimag = nreal*right.imaginary() + nimag*right.real();
}
public void divide(complex right )
{
double a=nreal*nreal+nimag*nimag;
nreal = ( nreal*right.real() + nimag*right.imaginary())/a;
nimag = (-nreal*right.imaginary() + nimag*right.real())/a;
}
public static complex add(complex left, complex right)
{ // return sum of two complex numbers
double r1=(left.real() + right.real());
double i1=(left.imaginary() + right.imaginary());
complex result;
result=new complex(r1,i1);
return result;
}
public static complex substract(complex left, complex right)
{ // return substraction of two complex numbers
complex result;
result=new complex((left.real() - right.real()),
(left.imaginary() - right.imaginary()));
return result;
}
public static complex multiply(complex left, complex right)
{ // return multiplication of two complex numbers
complex result;
result=new complex
((left.real()*right.real() - left.imaginary()*right.imaginary()),
(left.real()*right.imaginary() + left.imaginary()*right.real()));
return result;
}
public static complex divide(complex left, complex right)
{ // return division of two complex numbers
double a=right.real()*right.real()+right.imaginary()*right.imaginary();
complex result;
result=new complex
((left.real()*right.real() + left.imaginary()*right.imaginary())/a,
(-left.real()*right.imaginary() + left.imaginary()*right.real())/a);
return result;
}
public static complex pow(complex left, double right)
{ // return sum of two complex numbers
double Rad,th;
Rad=Math.pow(left.R(),right);
th=right*left.theta();
complex result;
result =new complex((Rad*Math.cos(th) ),
(Rad*Math.sin(th) ) );
return result;
}
public boolean smaller(complex left,complex right)
{
// less then comparison of two complex numbers
return (left.R() < right.R());
}
public boolean smaller_equal(complex left,complex right)
{
// less then and equal comparison of two complex numbers
return (left.R() <= right.R());
}
public boolean greater(complex left,complex right)
{
// greater then comparison of two complex numbers
return left.R() > right.R();
}
public boolean greater_equal(complex left,complex right)
{
// greater then and equal comparison of two complex numbers
return left.R() >= right.R();
}
public boolean equal(complex left,complex right)
{
// equal comparison of two complex numbers
return left.R() == right.R();
}
public boolean not_equal(complex left,complex right)
{
// not equal comparison of two complex numbers
return left.R() != right.R();
}
public static String toString(complex value)
{
String b="";
if(Math.abs(value.imaginary())!=1)
{
if(value.imaginary() >= 0)
b=b+"("+value.real()+" + "+value.imaginary()+"i )";
else
b=b+"("+value.real()+" - "+(-value.imaginary())+"i )";
}
else
{
if(value.imaginary() >= 0)
b=b+"("+value.real()+" + i )";
else
b=b+"("+value.real()+" - i )";
}
return b;
}
public String toString()
{
String b="";
if(Math.abs(imaginary())!=1)
{
if(imaginary() > 0)
b=b+"("+real()+" + "+imaginary()+"i )";
else if(imaginary() <0 )
b=b+"("+real()+" - "+(-imaginary())+"i )";
else
b=b+real()+ " ";
}
else
{
if(imaginary() > 0)
b=b+"("+real()+" + i )";
else if(imaginary() < 0)
b=b+"("+real()+" - i )";
else
b=b+real()+" ";
}
return b;
}
// data areas
double nreal;
double nimag;
};
//end of class complex definations
Şimdi Matrix metoduna daha detaylı bakabiliriz. İlkönce Matrix sınıfının metot başlıklarını verelim :
toString Metotları
public static String toString(double left)
left double değişkeninin String eşdeğerini yazıma hazır hale getirir.
public static String toString(double[][] left)
left double matrisini
String değişkeni olarak yazıma hazır hale getirir.
public static String toString(complex[][] left)
left complex
matrisini String değişkeni olarak yazıma hazır hale getirir.
public static String toStringT(complex[] left)
left complex
vektörünü String değişkeni sütun olarak yazıma hazır hale getirir.
public static String toString(complex[] left)
left complex
vektörünü String değişkeni satır olarak yazıma hazır hale getirir.
public static String toStringT(double[] left)
left double vektörünü
String değişkeni sütun olarak yazıma hazır hale getirir.
public static String toString(double[] left)
left complex
vektörünü String değişkeni satır olarak yazıma hazır hale getirir.
inverse matris ve denklem sistemi çözüm hesapları
Denklem çözümlerinde temel olarak iki sayısal analiz tekniğinden yararlanılmıştır. Birinci metod tam pivotlu gauss eliminasyon yöntem,, diğeri ise LU ayrıştırma yöntemidir. Gauss eliminasyon yöntemi matrisi eliminasyon yöntemi kullanarak bir üst üçgen matris haline dönüştürür, yani diagonal elementlerin altında kalan elemanları sıfırlar. Böylece denklem sistemi tek bilinmiyenli bir lineer denklemden başayarak iki bilinmiyenli, üç bilinmiyenli vs bir denklem sistemine dönüşür. Bu denklem sistemi bir bilinmiyenli denklemin çözülüp değerinin iki blinmiyenli denklemde yerine konulmasıyla ikincinin çözülmesi ve aynı işlemin seri halde devamıyla tüm bilinmeyenler çözülür. LU ayrıştırma metodunda ise matris önce bir üst üçgen matris(U) ve bir alt üçgen matris (L) ye ayrıştırılır. Sistem L ve U matrislerini kullanarak çözülür. İnverse (ters) matris hesabı da yine bu iki yöntemin kullanılmasıyla gerçekleştirilir. Eğer sadece bir matris çözülecekse gauss eliminasyon ve LU yöntemleri arasında hesaplama verimi farkı yoktur, ancak sadece denklemin sağ tarafı değişiyorsa LU daha verimli olarak hesap yapar. Hesaplama hassasiyetini arttırmak için pivoting adıverilen bir metod kullanılmıştır. Pivoting matrisin satır ve sütunlarının yerini değiştirerek yuvarlama hatalarının minimize edilmesi tekniğidir.
public static double[][] inv(double[][] a)
inverse
matrix hesaplar (tam pivotlu Gauss eliminasyon yöntemi)
public static double[][] inverse(double a[][])
inverse matrix hesaplar (tam pivotlu Gauss eliminasyon yöntemi)
public static double[] LUaxb(double a[][],double x[],int indx[])
lineer denklem sistemini LU metodu yardımıyla çözer.
public static double[] AXB(double a[][],double b[])
lineer denklem sistemini
tam pivotlu gauss eliminasyon metodu yardımıyla çözer.
Determinant hesapları
Determinant özel bir matris çarpım prosesidir. Tek bir sayıyla ifade edilir. Determinant hesabı için LU metodu daha uygun olduğundan bu metod kullanılmıştır. Aşağıdaki metodlar aslında tek bir metod olup, değişik isimleri verebilmek amacıyla üç ayrı metod olarak sunulmuştur.
public static double det(double a[][])
LU metodu
yardımıyla determinant hesaplar
public static double determinant(double a[][])
LU metodu
yardımıyla determinant hesaplar
public static double D(double a[][])
LU metodu
yardımıyla determinant hesaplar
Matrislerin dört işlemi
Burada matrislerin çarpımı (multiply) metodlarında vektörel matris çarpımlarından bahsediyoruz.
public static double[][]
multiply(double[][] left,double[][] right)
iki double matris çarpımı
public static double[] multiply(double[][] left,double[] right)
double matris - vektör çarpımı
public static double[] multiply(double[] left,double[][] right)
double
vektör - matris çarpımı
public static double[][] multiply(double left,double[][] right)
double
sayı - matris çarpımı
public static double[][] multiply(double[][] left,double right)
double
matris - sayı çarpımı
public static double[] multiply(double left,double[] right)
double
sayı - vektör çarpımı
public static double[] multiply(double[] left,double right)
double
sayı - vektör çarpımı
Matrislerin toplamı(add) aynı konumda olan elemanların toplamıdır. Fark
işlemi de (subtract) aynı konumda olan iki matrisden ilk verilenden ikinci
verilenin işaratinin değiştirilerek toplanmasıdır.
public static double[][] add(double[][] left,double[][] right)
iki double
matris toplamı
public static double[] add(double[] left,double[] right)
iki double
vektör toplamı
public static double[][] substract(double[][] left,double[][] right)
iki double
mtrisin farkı
public static double[] substract(double[] left,double[] right)
iki double
vektörün farkı
vectorlerin bölümüyle birinci vectorün tersinin (inverse) alınarak ikinci
matrisle çarpılması anlaşılır.
A/B = B*inv(A)
public static double[][] divide(double[][] left,double[][] right)
iki matrisin
bölümü (birinci matrisin ikinci matrisin inversüyle çarpımı)
public static double[][] LUdivide(double[][] left,double[][] right)
iki matrisin
LU metodu kullanılarak bölümü (birinci matrisin ikinci matrisin inversüyle
çarpımı)
public static double[] LUdivide(double[] left,double[][] right)
vektörün
matrise bölümü (LU metodu)
public static double[] divide(double[] left,double[][] right)
vektörün matrise bölümü
public static double[] divide(double[] left,double[][] right)
vektörün matrise bölümü
Üst işlemi : üst işleminden anlaşılan matristeki her elemanın ayrı ayrı üssünü almaktır
public static double[][] pow(double[][] right,double left)
matrisin double
kuvveti (her elemanının ayrı ayrı kuvveti)
public static double[] pow(double[] right,double left)
vektörün double kuvveti ( (her elemanının ayrı ayrı kuvveti)
Absulute matris : absolute matris, matris elemanlarının tümünün değerini positife çevirir.
public static double abs(double[][] left)
matrisin
mutlak değeri
public static double abs(double[] left)
vektörün mutlak değeri
Transpose matris : transpose matris matrisin satır ve sütünlarının yerini değiştirir, yani satırlar sütun, sütunlar satır haline dönüşür. (aij) T =aji
public static double[][] Transpose(double [][] left)
Transpose
matrix (satır ve sütunların yer değiştirmiş hali)
public static double[][] T(double [][] left)
Transpose
matrix (satır ve sütunların yer değiştirmiş hali)
public static complex[][] T(complex [][] left)
Transpose matrix (satır ve sütunların yer değiştirmiş hali)
Birim matris : birim matris diagonal elemanında (aii) 1 değeri olan diğer elemanları ise 0 taşıyan (aij,j!=I) matristir.
public static double[][] I(int n)
Birim matris
public static double[] one(int n)
Birim vektör (tüm değerler 1'e eşit)
EigenValue Problemi : EigenValue, veya sınır değer problemi matris işlemlerindeki ve matematikteki oldukça önemli bir işlemdir. Sınır değer işlemine kısaca bir denklem sisteminn köklerini bulma da diyebiliriz. Genel tanım olarak [A] matrisinin sınır değeri
[A]{X}- l[I]=0 olarak tanımlanır. Burada {X} denklemi sağlayan vektör setidir, eigen vektörü (sınır değer vektörü) olarak adlandırılır. l değeri ise EigenValue(sınır değeri) adını alır. Buradaki [I] birim matrisdir. Şüphesiz, verilen bir matris için birden fazla eigenvektör ve eigenvalue bulunabilir.
public static double[][]
balance(double b[][])
balance işlemi eigenvalue hesaplamalarını kolaylaştıran bir
ilk iterasyon prosesidir. Bazı durumlarda sonuca
ulaşmayabilir. Bu yüzden dikkatli kullanılmalıdır.
public static double[][] Hessenberg(double b[][])
Hessenberg
matris dönüşümü, eigenvalue hesaplamalarında kullanılır.
public static double[][] QR(double b[][])
Eigenvalue (sınır-değer hesaplama metodu, simetrik olmayan matrislere de uygulanabilir.)
sonuçlar n*2
(gerçek ve imajineri) matris olarak aktarılır.
public static double[][] eigenValue(double b[][])
QR metodunun
aynıdır.
public static complex[] eigenValueC(double b[][])
QR metodunun
aynıdır tek farkı sonuçları kompleks vektör olarak vermesidir.
public static double[][] balancedEigenValue(double b[][])
QR metodu
uygulanmadan önce balance metodunu uygular.
public static double[][] eigenQL(double a[][])
Eigenvalue
(sınır-değer hesaplama metodu, sadece simetrik matrislere uygulanabilir,
simetrik eigenvalue değerleri gerçektir,, komplex kökler vermez)
Polinomun kökleri
Polinomun köklerinin bulunması aslında bir eigenvalue (sınırdeğer) problemidir. Herhangi bir dereceden polinomun katsayıları verildiğinde karekteristikmatris adı verilen bir matris eşdeğeri oluşturulabilir. Bu matrisin eigenvalue değerleri verilen polinomun kökleridir.
public static double[][] poly_roots(double c[])
QR Eigenvalue (sınır-değer hesaplama metodu, simetrik olmayan matrislere de uygulanabilir.)
metodu kullanarak polinomun
köklerini hesaplar. polinom katsayılarını önce karektteristik matrise çevirir.
sonuçlar n*2 (gerçek ve imajineri) matris olarak aktarılır.
public static complex[] poly_rootsC(double c[])
Yukarıdakinin aynıdır. ancak sonuçlar kompleks vektör olarak aktarılır.
12.1.1 Matrix Sınıfının Listesi :
Problem 12.2 Matrix.java
, Matris proseslerini tanımlayan Matrix sınıfı
//=============================================================
// Numerical Analysis Packages in Java
// Matrix calculations class Matrix
// Author : Dr. Turhan Coban
//=============================================================
//Turhan Coban
import java.io.*;
import complex;
public class Matrix
{
// This class defines matrix and vector
// calculation methods
public static String
toString(double left)
{
//arrange double to string conversion so that all the
//matrix double variables nicely printed in the same column
String s="";
if(left=0) s=s+" ";
s=s+left;
double n=s.length();
while(n<13)
{
s=s+" ";
n=s.length();
}
return s;
}
public static String
toString(double[][] left)
{
//return a string representation of a matrix
int n,m;
String b;
b="";
n=left.length;
m=left[0].length;
for(int i=0;i<n;i++)
{
for(int j=0;j<m;j++)
{
b=b+toString(left[i][j]);
}
b=b+"\n";
}
return b;
}
public static String
toString(complex[][] left)
{
//return a string representation of a complex matrix
int n,m;
String b;
b="";
n=left.length;
m=left[0].length;
for(int i=0;i<n;i++)
{
for(int j=0;j<m;j++)
{
b=b+left[i][j].toString()+"\t";
}
b=b+"\n";
}
return b;
}
public static String
toStringT(complex[] left)
{
// returns a horizontal string representation of
// a complex vector
int n,m;
String b;
b="";
n=left.length;
for(int i=0;i<n;i++)
{
b=b+left[i].toString()+"\n";
}
return b;
}
public static String
toString(complex[] left)
{
// returns a vertical string representation of
// a complex vector
int n,m;
String b;
b="";
n=left.length;
for(int i=0;i<n;i++)
{
b=b+left[i].toString()+" ";
}
b=b+"\n";
return b;
}
public static String
toStringT(double[] left)
{
// returns a vertical string representation
// of a double vector
int n,m;
String b;
b="";
n=left.length;
for(int i=0;i<n;i++)
{
b=b+toString(left[i])+"\n";
}
return b;
}
public static String
toString(double[] left)
{
// returns a horizontal string representation
// of a double vector
int n,m;
String b;
b="";
n=left.length;
for(int i=0;i<n;i++)
{
b=b+toString(left[i]);
}
b=b+"\n";
return b;
}
public static String toString(int[] left)
{
// returns a horizontal string representation
// of a integer vector
int n,m;
String b;
b="";
n=left.length;
for(int i=0;i<n;i++)
{
b=b+left[i]+"\t";
}
b=b+"\n";
return b;
}
public static double
SIGN(double a,double b)
{
//returns the value of double a with sign of double b;
//if a=-2, b= 3 SIGN(a,b) returns 2
//if a=-2, b=-3 SIGN(a,b) returns -2
//if a= 2, b=-3 SIGN(a,b) returns -2
//if a= 2, b= 3 SIGN(a,b) returns 2
if(b!=0)
return Math.abs(a)*b/Math.abs(b);
else
return Math.abs(a);
}
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[][]
inverse(double a[][])
{
//inverse of a matrix
//this method enable usage of inv or inverse as
//name of mtrix inversion method
return Matrix.inv(a);
}
//LU decomposition method
public static double[][] LU(double c[][],int indx[],int d[])
{
//returns LU decomposition of matrix c and index indx
double a[][];
int n=c.length;
a=new double[n][n];
double vv[];
vv=new double[n];
double sum,dum,big,temp;
int i,j,k;
int imax;
int nmax=100;
double tiny=1.0e-40;
imax=0;
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
a[i-1][j-1]=c[i-1][j-1];
}
d[0]=1;
for(i=1;i<=n;i++)
{
big=0.0;
for(j=1;j<=n;j++)
{
if(Math.abs(a[i-1][j-1])big) big=Math.abs(a[i-1][j-1]);
}
if(big==0) {System.out.println("singular matrix");return a;}
vv[i-1]=1.0/big;
}
for(j=1;j<=n;j++)
{
for(i=1;i<j;i++)
{
sum=a[i-1][j-1];
for(k=1;k<i;k++)
{
sum-=a[i-1][k-1]*a[k-1][j-1];
}
a[i-1][j-1]=sum;
}
big=0;
for(i=j;i<=n;i++)
{
sum=a[i-1][j-1];
for(k=1;k<j;k++)
{
sum-=a[i-1][k-1]*a[k-1][j-1];
}
a[i-1][j-1]=sum;
dum=vv[i-1]*Math.abs(sum);
if(dum=big)
{
imax=i;
big=dum;
}
} //end of i=0
if(j != imax)
{
for(k=1;k<=n;k++)
{
dum=a[imax-1][k-1];
a[imax-1][k-1]=a[j-1][k-1];
a[j-1][k-1]=dum;
}
d[0]=-d[0];
vv[imax-1]=vv[j-1];
} //end of if
indx[j-1]=imax;
if(a[j-1][j-1]==0) a[j-1][j-1]=tiny;
if(j!=n)
{
dum=1.0/a[j-1][j-1];
for(i=j+1;i<=n;i++)
a[i-1][j-1]*=dum;
}//endif
} //end for j=
return a;
}
public static double[]
LUaxb(double a[][],double x[],int indx[])
{
//solves AX=B system of linear equation of LU decomposed matrix a
//(calculated by method LU)
int ii=0;
int i,j,ll=0;
double sum=0;
int n=a.length;
double b[];
b=new double[n];
for(i=1;i<=n;i++)
{
b[i-1]=x[i-1];
}
for(i=1;i<=n;i++)
{
ll=indx[i-1];
sum=b[ll-1];
b[ll-1]=b[i-1];
if(ii!=0)
{
for(j=ii;j<=(i-1);j++)
{
sum-=a[i-1][j-1]*b[j-1];
}
}
else if(sum!=0) ii=i;
b[i-1]=sum;
}
for(i=n;i=1;i--)
{
sum=b[i-1];
if(i<n)
{
for(j=(i+1);j<=n;j++)
{
sum-=a[i-1][j-1]*b[j-1];
}
}
b[i-1]=sum/a[i-1][i-1];
}
return b;
}
public static double[]
AXB(double a[][],double b[])
{
//Solution of system of linear equations by LU method
// note that the same calculation can be done by divide method.
int n=a.length;
double c[]=new double[n];
int d[]={1};
int indx[]=new int[n];
double e[][]=new double[n][n];
e=Matrix.LU(a,indx,d);
c=Matrix.LUaxb(e,b,indx);
return c;
} //end of AXB
public static double[][]
LUinv(double a[][])
{
//inverse of a matrix by using LU decomposition method
//this method is more efficient than inv (or inverse)
int n=a.length;
double c[][]=new double[n][n];
double b[][]=Matrix.I(n);
int d[]={0};
int indx[]=new int[n];
double e[][]=new double[n][n];
e=Matrix.LU(a,indx,d);
for(int i=0;i<n;i++)
{
c[i]=Matrix.LUaxb(e,b[i],indx);
}
return Matrix.T(c);
} //end of LUinv
public static double
det(double a[][])
{
//determinant of a matrix
int n=a.length;
int indx[]=new int[n];
int d[]={1};
double e;
double b[][]=new double[n][n];
b=Matrix.LU(a,indx,d);
e=d[0];
for(int i=0;i<n;i++)
e=e*b[i][i];
return e;
} //end of det
public static double
determinant(double a[][])
{
//determinant of a matrix
return Matrix.det(a);
}
public static double
D(double a[][])
{
//determinant of a matrix
return Matrix.det(a);
}
//************* multiply
methods definitions *********************
public static double[][] multiply(double[][] left,double[][] right)
{
//multiplication of two matrices
int ii,jj,i,j,k;
int m1=left[0].length;
int n1=left.length;
int m2=right[0].length;
int n2=right.length;
double[][] b;
b=new double[m1][n2];
if(n1 != m2)
{
System.out.println("inner matrix dimensions must agree");
for(ii=0;ii<n1;ii++)
{
for(jj=0;jj<m2;jj++)
b[ii][jj]=0;
}
return b;
}
for(i=0;i<m1;i++)
{
for(j=0;j<n2;j++)
{
for(k=0;k<n1;k++)
b[i][j]+=left[i][k]*right[k][j];
}
}
return b;
//end of multiply of two matrices
}
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[]
multiply(double[] left,double[][] right)
{
//multiplication of one vector with one matrix
int ii,jj,i,j,k;
int m2=right[0].length;
int n2=right.length;
int m1=left.length;
double[] b;
b=new double[m1];
if(n2 != m1)
{
System.out.println("inner matrix dimensions must agree");
for(ii=0;ii<n2;ii++)
{
b[ii]=0;
}
return b;
}
for(i=0;i<m2;i++)
{
b[i]=0;
for(k=0;k<m1;k++)
b[i]+=right[i][k]*left[k];
}
return b;
//end of multiply of a vector and a matrix
}
public static double[][]
multiply(double left,double[][] right)
{
//multiplying a matrix with a constant
int i,j;
int n=right.length;
int m=right[0].length;
double b[][];
b=new double[n][m];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
b[i][j]=right[i][j]*left;
}
return b;
//end of multiplying a matrix with a constant double
}
public static double[][]
multiply(double[][] left,double right)
{
//multiplying a matrix with a constant
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[n][m];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
b[i][j]=left[i][j]*right;
}
return b;
//end of multiplying a matrix with a constant double
}
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)
{
//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;
}
//*************** end of multiply methods definitions **************
//=============== defination of power methods pow ==============
public static double[][] pow(double[][] right,double left)
{
// power of a matrix
int i,j;
double b[][];
int n=right.length;
int m=right[0].length;
b=new double[n][m];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
{
if(left==0.0)
{
b[i][j]=1.0;
}
else
{
b[i][j]=Math.pow(right[i][j],left);
}
}
}
return b;
//end of power of a matrix
}
public static double[]
pow(double[] right,double left)
{
// power of a vector
int i;
int n=right.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
{
if(left==0.0)
{
b[i]=1.0;
}
else
{
b[i]=Math.pow(right[i],left);
}
}
return b;
//end of power of a vector
}
//=================end
of power method pow definitions =============
//***************** addition add methods **************************
public static double[][]
add(double[][] left,double[][] right)
{
//addition of two matrices
int n1=left.length;
int m1=left[0].length;
int n2=right.length;
int m2=right[0].length;
int nMax,mMax;
int i,j;
if(m1=m2) mMax=m1;
else mMax=m2;
if(n1=n2) nMax=n1;
else nMax=n2;
double b[][];
b=new double[nMax][mMax];
for(i=0;i<n1;i++)
{
for(j=0;j<m1;j++)
{
b[i][j]=b[i][j]+left[i][j];
}
}
for(i=0;i<n2;i++)
{
for(j=0;j<m2;j++)
{
b[i][j]=b[i][j]+right[i][j];
}
}
return b;
//end of matrix addition method
}
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[][]
substract(double[][] left,double[][] right)
{
//addition of two matrices
int n1=left.length;
int m1=left[0].length;
int n2=right.length;
int m2=right[0].length;
int nMax,mMax;
int i,j;
if(m1=m2) mMax=m1;
else mMax=m2;
if(n1=n2) nMax=n1;
else nMax=n2;
double b[][];
b=new double[nMax][mMax];
for(i=0;i<n1;i++)
{
for(j=0;j<m1;j++)
{
b[i][j]=b[i][j]+left[i][j];
}
}
for(i=0;i<n2;i++)
{
for(j=0;j<m2;j++)
{
b[i][j]=b[i][j]-right[i][j];
}
}
return b;
//end of matrix substraction method
}
public static double[]
substract(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 substraction method
}
//============== division
of the matrices
public static double[][] divide(double[][] left,double[][] right)
{
//division of two matrices
int n=right.length;
int m=right[0].length;
double b[][];
b=new double[n][m];
b=Matrix.multiply(Matrix.inv(right),left);
return b;
}
public static double[][]
LUdivide(double[][] left,double[][] right)
{
//division of two matrices utilises LUinv method instead of inv
int n=right.length;
int m=right[0].length;
double b[][];
b=new double[n][m];
b=Matrix.multiply(Matrix.LUinv(right),left);
return b;
}
public static double[]
divide(double[] left,double[][] right)
{
//division of two matrices
int n=right.length;
int m=right[0].length;
double b[];
b=new double[n];
b=Matrix.multiply(Matrix.inv(right),left);
return b;
}
public static double[]
LUdivide(double[] left,double[][] right)
{
//division of two matrices utilises AXB (LU decomposition method)
//in fact this method is exactly same as AXB except spacing of the
//arguments
return AXB(right,left);
}
//============== absolute value of a matrix===================
public static double
abs(double[][] left)
{
// absoulute value of a matrix
int i,j;
int n=left.length;
int m=left[0].length;
double b=0;
for(i=0;i<n;i++)
for(j=0;j<m;j++)
{
b=b+Math.abs(left[i][j]);
}
return b;
}
public static double
abs(double[] left)
{
// absolute value of a vector
int i;
int n=left.length;
double b=0;
for(i=0;i<n;i++)
{
b=b+Math.abs(left[i]);
}
return b;
}
//===============special
matrices==============================
public static double[][] Transpose(double [][] left)
{
//transpose matrix (if A=a(i,j) Transpose(A)=a(j,i)
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[m][n];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
{
b[j][i]=left[i][j];
}
}
return b;
}
public static double[][]
T(double [][] left)
{
//transpose matrix (if A=a(i,j) T(A)=a(j,i)
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[m][n];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
{
b[j][i]=left[i][j];
}
}
return b;
}
public static complex[][]
T(complex [][] left)
{
//transpose matrix (if A=a(i,j) T(A)=a(j,i)
int i,j;
int n=left.length;
int m=left[0].length;
complex b[][];
b=new complex[m][n];
for(i=0;i<n;i++)
{
for(j=0;j<m;j++)
{
b[j][i]=new complex(left[i][j]);
}
}
return b;
}
public static double[][]
I(int n)
{
//unit matrix
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
b[i][i]=1.0;
return b;
}
public static double[]
one(int n)
{
//one matrix
double b[];
b=new double[n];
for(int i=0;i<n;i++)
b[i]=1.0;
return b;
}
public static double[][]
characteristic_matrix(double c[])
{
//this routine converts polynomial coefficients to a matrix
//with the same eigenvalues (roots)
int n=c.length-1;
int i;
double a[][]=new double[n][n];
for(i=0;i<n;i++)
{
a[0][i]=-c[i+1]/c[0];
}
for(i=0;i<n-1;i++)
{
a[i+1][i]=1;
}
return a;
}
//===========Eigen value calculations ==============
public static double[][]
balance(double b[][])
{
// balance of a matrix for more accurate eigenvalue
// calculations
double radix=2.0;
double sqrdx=radix*radix;
double c,r,f,s,g;
int m,j,i,last;
int n=b.length;
last=0;
double a[][];
a=new double[n][n];
f=1;
s=1;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
a[i-1][j-1]=b[i-1][j-1];
while(last==0)
{
last=1;
for(i=1;i<=n;i++)
{
c=0;r=0;
for(j=1;j<=n;j++)
{
if(j != i)
{
c+=Math.abs(a[j-1][i-1]);
r+=Math.abs(a[i-1][j-1]);
} //end of if(j!=..
} //end of for(j=1...
if(c != 0 && r != 0 )
{
g=r/radix;
f=1.0;
s=c+r;
while(c<g)
{
f*=radix;
c*=sqrdx;
}
g=r*radix;
while(cg)
{
f/=radix;
c/=sqrdx;
}
} //end of if(c != 0 && ....
if( (c+r)/f < 0.95*s )
{
last=0;
g=1.0/f;
for(j=1;j<=n;j++) { a[i-1][j-1]*=g; }
for(j=1;j<=n;j++) { a[j-1][i-1]*=f; }
} //end of if( ((c+r..
}//end of for(i=1;i<=n....
} //end of while last==0
return a;
}
public static double[][]
Hessenberg(double b[][])
{
// Calculates the hessenberg matrix
// it is used in QR method to calculate eigenvalues
// of a matrix(symmetric or non-symmetric)
int m,j,i;
int n=b.length;
double a[][];
a=new double[n][n];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i][j]=b[i][j];
double x,y;
if(n2)
{
for(m=2;m<=(n-1);m++)
{
x=0.0;
i=m;
for(j=m;j<=n;j++)
{
if(Math.abs(a[j-1][m-2]) Math.abs(x))
{
x=a[j-1][m-2];
i=j;
} //end of if(Math.abs(..
}//end of for(j=m,j<=n...
if(i!=m)
{
for(j=(m-1);j<=n;j++)
{
y=a[i-1][j-1];
a[i-1][j-1]=a[m-1][j-1];
a[m-1][j-1]=y;
}//end of for(j=(m-1)..
for(j=1;j<=n;j++)
{
y=a[j-1][i-1];
a[j-1][i-1]=a[j-1][m-1];
a[j-1][m-1]=y;
}//end of for(j=1;j<=n....
} //end of if(i!=m)
if(x != 0.0)
{
for(i=(m+1);i<=n;i++)
{
y=a[i-1][m-2];
if(y!=0.0)
{
y=y/x;
a[i-1][m-2]=y;
for(j=m;j<=n;j++)
{
a[i-1][j-1]-=y*a[m-1][j-1];
}
for(j=1;j<=n;j++)
{
a[j-1][m-1]+=y*a[j-1][i-1];
}
}//end of if(y!=0..
}//end of for(i=(m+1)...
} //end of if(x != 0.0...
}//end of for(m=2;m<=(n-1)..
}//end of Hessenberg
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
if(i(j+1)) a[i-1][j-1]=0;
}
return a;
}
public static double[][]
QR(double b[][])
{
//calculates eigenvalues of a Hessenberg matrix
int n=b.length;
double rm[][]=new double[2][n];
double a[][]=new double[n+1][n+1];
double wr[]=new double[n+1];
double wi[]=new double[n+1];
int nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r=0,q=0,p=0,anorm;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i+1][j+1]=b[i][j];
anorm=Math.abs(a[1][1]);
for (i=2;i<=n;i++)
for (j=(i-1);j<=n;j++)
anorm += Math.abs(a[i][j]);
nn=n;
t=0.0;
while (nn = 1) {
its=0;
do {
for (l=nn;l=2;l--) {
s=Math.abs(a[l-1][l-1])+Math.abs(a[l][l]);
if (s == 0.0) s=anorm;
if ((double)(Math.abs(a[l][l-1]) + s) == s) break;
}
x=a[nn][nn];
if (l == nn) {
wr[nn]=x+t;
wi[nn--]=0.0;
}
else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == (nn-1)) {
p=0.5*(y-x);
q=p*p+w;
z=Math.sqrt(Math.abs(q));
x += t;
if (q = 0.0) {
z=p+Matrix.SIGN(z,p);
wr[nn-1]=wr[nn]=x+z;
if (z!=0) wr[nn]=x-w/z;
wi[nn-1]=wi[nn]=0.0;
}
else {
wr[nn-1]=wr[nn]=x+p;
wi[nn-1]= -(wi[nn]=z);
}
nn -= 2;
}
else {
if (its == 30) System.out.println("Too many iterations in hqr");
if (its == 10 || its == 20) {
t += x;
for (i=1;i<=nn;i++) a[i][i] -= x;
s=Math.abs(a[nn][nn-1])+Math.abs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=(nn-2);m=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=Math.abs(p)+Math.abs(q)+Math.abs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=Math.abs(a[m][m-1])*(Math.abs(q)+Math.abs(r));
v=Math.abs(p)*(Math.abs(a[m-1][m-1])+
Math.abs(z)+Math.abs(a[m+1][m+1]));
if ((double)(u+v) == v) break;
}
for (i=m+2;i<=nn;i++) {
a[i][i-2]=0.0;
if (i != (m+2)) a[i][i-3]=0.0;
}
for (k=m;k<=nn-1;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k != (nn-1)) r=a[k+2][k-1];
if ((x=Math.abs(p)+Math.abs(q)+Math.abs(r)) != 0.0) {
p /= x;
q /= x;
r /= x;
}
}
if ((s=Matrix.SIGN(Math.sqrt(p*p+q*q+r*r),p)) != 0.0) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
}
else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<=nn;j++) {
p=a[k][j]+q*a[k+1][j];
if (k != (nn-1)) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn<k+3 ? nn : k+3;
for (i=l;i<=mmin;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k != (nn-1)) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
}
}
}
}
} while (l < nn-1);
}
for(i=0;i<n;i++)
{
rm[0][i]=wr[i+1];
rm[1][i]=wi[i+1];
}
return rm;
} //end of QR
public static double[][]
eigenValue(double b[][])
{
// this routine input a matrix (non symetric or symmetric)
// and calculate eigen values
// method balance can be used prior to this method to balance
// the input matrix
int n=b.length;
double d[][]=new double[2][n];
d=Matrix.QR(Matrix.Hessenberg(b));
return d;
}
public static complex[]
eigenValueC(double b[][])
{
// this routine input a matrix (non symetric or symmetric)
// and calculate eigen values
// method balance can be used prior to this method to balance
// the input matrix
//output eigenvalues will be in a vector of complex form
int n=b.length;
double d[][]=new double[2][n];
d=Matrix.QR(Matrix.Hessenberg(b));
complex c[]=new complex[n];
for(int i=0;i<n;i++)
{
c[i]=new complex(d[0][i],d[1][i]);
}
return c;
}
//roots of a polynomial
public static double[][] poly_roots(double c[])
{
//roots of a degree n polynomial
// P(x)=c[n]*x^n+c[n-1]*x^(n-1)+....+c[1]*x+c[0]=0;
int n=c.length-1;
double a[][]=new double[n][n];
a=characteristic_matrix(c);
double d[][]=new double[2][n];
d=balancedEigenValue(a);
return d;
}
public static complex[]
poly_rootsC(double c[])
{
// roots of a degree n polynomial
// P(x)=c[n]*x^n+c[n-1]*x^(n-1)+....+c[1]*x+c[0]=0;
// roots are returned as complex variables
int n=c.length-1;
double a[][]=new double[n][n];
a=characteristic_matrix(c);
double d[][]=new double[2][n];
d=balancedEigenValue(a);
complex e[]=new complex[n];
for(int i=0;i<n;i++)
e[i]=new complex(d[0][i],d[1][i]);
return e;
}
public static double[][]
balancedEigenValue(double b[][])
{
// this routine input a matrix (non symetric or symmetric)
// and calculates eigen values
// method balance is used to balance the matrix previous to
// actual calculations
int n=b.length;
double d[][]=new double[2][n];
d=Matrix.QR(Matrix.Hessenberg(Matrix.balance(b)));
return d;
}
public static double[][]
tridiagonal(double b[][], double d[], double e[])
{
//reduces matrix to tridiaonal form by using householder transformation
//this method is used by QL method to calculate eigen values
//and eigen vectors of a symmetric matrix
int l,k,j,i;
int n=b.length;
double scale,hh,h,g,f;
double a[][]=new double[n+1][n+1];
double c[][]=new double[n][n];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i][j]=b[i][j];
for (i=n;i=2;i--) {
l=i-1;
h=scale=0.0;
if (l 1) {
for (k=1;k<=l;k++)
scale += Math.abs(a[i-1][k-1]);
if (scale == 0.0)
e[i-1]=a[i-1][l-1];
else {
for (k=1;k<=l;k++) {
a[i-1][k-1] /= scale;
h += a[i-1][k-1]*a[i-1][k-1];
}
f=a[i-1][l-1];
g=(f = 0.0 ? -Math.sqrt(h) : Math.sqrt(h));
e[i-1]=scale*g;
h -= f*g;
a[i-1][l-1]=f-g;
f=0.0;
for (j=1;j<=l;j++) {
a[j-1][i-1]=a[i-1][j-1]/h;
g=0.0;
for (k=1;k<=j;k++)
g += a[j-1][k-1]*a[i-1][k-1];
for (k=j+1;k<=l;k++)
g += a[k-1][j-1]*a[i-1][k-1];
e[j-1]=g/h;
f += e[j-1]*a[i-1][j-1];
}
hh=f/(h+h);
for (j=1;j<=l;j++) {
f=a[i-1][j-1];
e[j-1]=g=e[j-1]-hh*f;
for (k=1;k<=j;k++)
a[j-1][k-1] -= (f*e[k-1]+g*a[i-1][k-1]);
}
}
} else
e[i-1]=a[i-1][l-1];
d[i-1]=h;
}
d[1-1]=0.0;
e[1-1]=0.0;
/* Contents of this loop can be omitted if eigenvectors not
wanted except for statement d[i-1]=a[i-1][i-1]; */
for (i=1;i<=n;i++) {
l=i-1;
if (d[i-1] != 0) {
for (j=1;j<=l;j++) {
g=0.0;
for (k=1;k<=l;k++)
g += a[i-1][k-1]*a[k-1][j-1];
for (k=1;k<=l;k++)
a[k-1][j-1] -= g*a[k-1][i-1];
}
}
d[i-1]=a[i-1][i-1];
a[i-1][i-1]=1.0;
for (j=1;j<=l;j++) a[j-1][i-1]=a[i-1][j-1]=0.0;
}
return a;
}
public static double
pythag(double a, double b)
{
//this method is used by QL method
double absa,absb;
absa=Math.abs(a);
absb=Math.abs(b);
if (absa absb) return absa*Math.sqrt(1.0+(absb/absa)*(absb/absa));
else return (absb==0.0 ? 0.0 : absb*Math.sqrt(1.0+(absa/absb)*(absa/absb)));
}
public static double[][]
QL(double d[], double e[], double a[][])
{
// QL algorithm : eigenvalues of a symmetric matrix reduced to tridiagonal
// form by using method tridiagonal
int n=d.length;
int m,l,iter,i,j,k;
double s,r,p,g,f,dd,c,b;
for (i=2;i<=n;i++) e[i-2]=e[i-1];
e[n-1]=0.0;
double z[][]=new double[n][n];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
z[i][j]=a[i][j];
for (l=1;l<=n;l++) {
iter=0;
do {
for (m=l;m<=n-1;m++) {
dd=Math.abs(d[m-1])+Math.abs(d[m]);
if ((double)(Math.abs(e[m-1])+dd) == dd) break;
}
if (m != l) {
if (iter++ == 30) System.out.println("Too many iterations
in QL");
g=(d[l]-d[l-1])/(2.0*e[l-1]);
r=Matrix.pythag(g,1.0);
g=d[m-1]-d[l-1]+e[l-1]/(g+Matrix.SIGN(r,g));
s=c=1.0;
p=0.0;
for (i=m-1;i=l;i--) {
f=s*e[i-1];
b=c*e[i-1];
e[i]=(r=Matrix.pythag(f,g));
if (r == 0.0) {
d[i] -= p;
e[m-1]=0.0;
break;
}
s=f/r;
c=g/r;
g=d[i ]-p;
r=(d[i-1]-g)*s+2.0*c*b;
d[i ]=g+(p=s*r);
g=c*r-b;
for (k=1;k<=n;k++) {
f=z[k-1][i ];
z[k-1][i ]=s*z[k-1][i-1]+c*f;
z[k-1][i-1]=c*z[k-1][i-1]-s*f;
}
}
if (r == 0.0 && i = l) continue;
d[l-1] -= p;
e[l-1]=g;
e[m-1]=0.0;
}
} while (m != l);
}
return z;
}
public static double[][]
eigenQL(double a[][])
{
// QL algoritm to solve eigen value problems
// symmetric matrices only (real eigen values)
// first column of the matrix returns eigen values
// second..n+1 column returns eigen vectors.
// Note : If matrix is not symmetric DO NOT use
// this method use eigenValue method (a QR algorithm)
int i,j;
int n=a.length;
double sum[]=new double[n];;
double d[]=new double[n];
double b[][]=new double[n][n];
double e[]=new double[n];
double z[][]=new double[n+1][n];
b=tridiagonal(a,d,e);
b=QL(d,e,b);
for(j=0;j<n;j++)
{
z[0][j]=d[j];
for(i=0;i<n;i++)
{
z[i+1][j]=b[i][j]/b[0][j];
if(z[i+1][j]<1e-13) z[i+1][i]=0;
}
}
return z;
}
//end of eigen value programs
//end of class Matrix
}
12.3 MATRIX SINIFI ÖRNEK PROGRAMLARI
İnverse
matris : Inverse
matrisi hesaplar. Buradaki temel hesaplama yolu tam pivotlu gauss eliminasyon
metodudur. Bu metod göreceli olarak küçük matrisler içindir. dev boyutlardaki
matrisler için tavsiye edilmez. Çok büyük boyutlarda iterativ metodlar
daha az sayıda işlemle çözüme gidebilir. Bu metod hassas sonuç verebilme
özelliğinden dolayı seçilmiştir. aynı zamanda LU (alt ve üst üçgen matris)
parçalama metoduyla da çözebiliriz. Bu metodun avantajı lineer bir denklem
sistemi çözerken eğer ikinci taraf sürekli değişiyorsa ortaya çıkar. tek
bir çözüm için ilave bir avantaj getirmez.
Problem 12.3 Matrix3.java , inverse matris problemi örneği
import java.io.*;
import Numeric;
class Matrix3
{
public static void main(String args[]) throws IOException
{
/*
double a[][];
a=new double[5][5];
a[0][0]=1;
a[0][1]=2;
a[0][2]=0;
a[0][3]=0;
a[0][4]=0;
a[1][0]=-2;
a[1][1]=3;
a[1][2]=0;
a[1][3]=0;
a[1][4]=0;
a[2][0]=3;
a[2][1]=4;
a[2][2]=50;
a[2][3]=0;
a[2][4]=0;
a[3][0]=-4;
a[3][1]=5;
a[3][2]=-60;
a[3][3]=7;
a[3][4]=0;
a[4][0]=-5;
a[4][1]=6;
a[4][2]=-70;
a[4][3]=8;
a[4][4]=-9;
*/
double c[][]=new double[5][5];
double a[][]={{1,2,0,0,0},{-2,3,0,0,0},{3,4,50,0,0},{-4,5,-60,7,0},
{-5,6,-70,8,-9}};
double b[]=new double[5];
int d[]=new int[1];
b[0]=1;
b[1]=0;
b[2]=0;
b[3]=0;
b[4]=0;
System.out.println("Original Matrix : ");
System.out.println(Matrix.toString(a));
System.out.println("Inverse Matrix : (Method inv) ");
System.out.println(Matrix.toString(Matrix.inv(a)));
System.out.println("Matrix * Inverse Matrix : (Method multiply) ");
System.out.println(Matrix.toString(Matrix.multiply(a,Matrix.inv(a))));
System.out.println("Solution of system of equation : ");
System.out.println("with second side : ");
System.out.println(Matrix.toStringT(b));
System.out.println(Matrix.toString(Matrix.divide(b,a)));
int indx[]=new int[5];
d[0]=0;
c=Matrix.LU(a,indx,d);
System.out.println("LU decomposed matrix : (Method LU) ");
System.out.println(Matrix.toString(c));
System.out.println("Matrix inversion by LU decomposition : (Method
LUinv)");
System.out.println(Matrix.toString(Matrix.LUinv(a)));
System.out.println("Matrix * Inverse Matrix (LU Decomposition) :");
System.out.println(Matrix.toString(Matrix.multiply(a,Matrix.LUinv(a))));
}
}
inverse matrix metotları sonucu :
Original
Matrix :
1.0 2.0 0.0
0.0 0.0
-2.0 3.0 0.0 0.0 0.0
3.0 4.0 50.0 0.0 0.0
-4.0 5.0 -60.0 7.0 0.0
-5.0 6.0 -70.0 8.0 -9.0
Inverse Matrix : (Method inv)
0.4285714285714286 -0.2857142857142857 0.0 -1.1895246692412392E-17
3.96508223080413E-18
0.2857142857142857 0.14285714285714285 0.0 5.947623346206196E-18-1.9825411154020644E-18
-0.048571428571428564 0.005714285714285717 0.019999999999999997-3.8659551750340293E-19
1.898283117997477E-18
-0.37551020408163255-0.21632653061224483 0.1714285714285714 0.14285714285714282
2.2431036619977644E-17
-0.003628117913832162 0.01723356009070296-0.0031746031746031824 0.12698412698412698-0.11111111111111112
Matrix * Inverse Matrix : (Method multiply)
1.0 0.0 0.0 0.0 1.5407439555097887E-33
-1.1102230246251565E-16 1.0 0.0 4.163336342344337E-17-1.3877787807814454E-17
4.440892098500626E-16 1.6653345369377348E-16 0.9999999999999999-3.122502256758254E-17
9.887923813067798E-17
0.0 2.220446049250313E-16 0.0 0.9999999999999999 1.7347234759768022E-17
-2.8449465006019636E-16 1.942890293094024E-16 1.734723475976807E-16 0.0
1.0
Solution of system of equation :
with second side :
1.0
0.0
0.0
0.0
0.0
0.4285714285714286
0.2857142857142857
-0.048571428571428564
-0.37551020408163255
-0.003628117913832162
LU decomposed matrix : (Method LU)
-2.0 3.0 0.0 0.0 0.0
-0.5 3.5 0.0 0.0 0.0
2.5 -0.42857142857142855-70.0 8.0 -9.0
-1.5 2.4285714285714284-0.7142857142857143 5.714285714285714-6.428571428571429
2.0 -0.2857142857142857 0.8571428571428571 0.025000000000000064 7.875
Matrix inversion by LU decomposition : (Method LUinv)
0.42857142857142855-0.2857142857142857 -0.0 -0.0 -0.0
0.2857142857142857 0.14285714285714285 0.0 0.0 0.0
-0.04857142857142856 0.005714285714285712 0.02 -0.0 -0.0
-0.37551020408163255-0.21632653061224494 0.17142857142857143 0.14285714285714285
0.0
-0.0036281179138321806 0.017233560090702916-0.003174603174603183 0.12698412698412698-0.1111111111111111
Matrix * Inverse Matrix (LU Decomposition) :
1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0
4.440892098500626E-16-1.1102230246251565E-16 1.0 0.0 0.0
-4.440892098500626E-16-4.440892098500626E-16 0.0 1.0 0.0
-1.1796119636642288E-16-8.326672684688674E-17-4.5102810375396984E-17
0.0 1.0
Problem 12.4 Matrix5.java, inverse matris problemi örneği
import java.io.*;
import Numeric;
class Matrix5
{
public static void main(String args[]) throws IOException
{
//calculation of an inverse matrix
double a[][];
a=new double[2][2];
a[0][0]=1;
a[0][1]=2;
a[1][0]=3;
a[1][1]=4;
System.out.println("Matrix : ");
System.out.println(Matrix.toString(a));
System.out.println("Inverse Matrix : ");
System.out.println(Matrix.toString(Matrix.inv(a)));
System.out.println("Matrix * Inverse Matrix :");
System.out.println(Matrix.toString(Matrix.multiply(a,Matrix.inv(a))));
System.out.println("Determinant = "+Matrix.det(a));
}
}
Matrix :
1.0 2.0
3.0 4.0
Inverse Matrix :
-2.0 1.0
1.5 -0.5
Matrix * Inverse Matrix :
1.0 0.0
0.0 1.0
Determinant = -2.0
EigenValue - Sınır değer problemi : matislerle ilgili belki de en ilginç proses eigen value - sınır değer problemidir. Sınır değer hesaplamaları iteratif bir problemdir. Simetrik matrislerde sonuçlar gerçek sayı olarak çıkar. bu tür uygulamalarda QL metodu kullanılabilir. Simetrik olmayan matrislerin sınır değerleri kompleks değerlerdir (kompleks ve gerçek değerlerin karışımı olabilir.) burada temel olarak QR formülü kullanılmıştır. Eigen value problemleriyle direk olarak ilgili ilginç bir problem de ninci dereceden bir polinomun köklerinin hesaplanması problemidir. polinom karakteristik matris ismi verilen özel bir matrise dönüştürülebilir. Bu matrisin eigen-value değerleri polinomun da kökleridir.
Karekterisitik matrisin tanımı
public static double[][]
characteristic_matrix(double c[])
{
int n=c.length-1; int i;
double a[][]=new double[n][n];
for(i=0;i<n;i++)
{ a[0][i]=-c[i+1]/c[0]; }
for(i=0;i<n-1;i++)
{ a[i+1][i]=1; }
return a; }
metodu ile verilmiştir.
Problem 12.5 Matrix1.java , eigenvalue problemi
import java.io.*;
import Numeric;
class Matrix1
{
public static void main(String args[]) throws IOException
{
//Eigen Value calculations
double c[][];
c=new double[5][5];
double a[][]={{1,2,0,0,0}, {-2,3,0,0,0}, {3,4,50,0,0},
{-4,5,-60,7,0},{-5,6,-70,8,-9}};
System.out.println("Matris : ");
System.out.println(Matrix.toString(a));
complex c1[];
int n=a.length;
c1=new complex[n];
c1=Matrix.eigenValueC(a);
System.out.println("Kompleks Eigenvalue vektörü : ");
System.out.println(Matrix.toStringT(c1));
}
}
Matris :
1.0 2.0 0.0 0.0 0.0
-2.0 3.0 0.0 0.0 0.0
3.0 4.0 50.0 0.0 0.0
-4.0 5.0 -60.0 7.0 0.0
-5.0 6.0 -70.0 8.0 -9.0
Kompleks Eigenvalue vectörü :
50.0
(2.000000000000002 - 1.732050807568873i )
(2.000000000000002 + 1.732050807568873i )
7.000000000000001
-8.999999999999998
Problem 12.6 Matrix4.java , ninci dereceden denklemin kökleri ve eigenvalue problemleri
import java.io.*;
import Numeric;
class Matrix4
{
public static void main(String args[]) throws IOException
{
//Eigen Value calculations using QL method
double p[]={-1, 8, -19, 12};
double a[][]=Matrix.characteristic_matrix(p);
System.out.println("Orijinal Polinom : ");
System.out.println(Matrix.toString(p));
System.out.println("Karekteristik matris : ");
System.out.println(Matrix.toString(a));
System.out.println("Eigen Value (sinir deger) : ");
System.out.println(Matrix.toString(Matrix.eigenValue(a)));
System.out.println("Eigen Value (sinir deger) balancedeigenValue metodu : ");
System.out.println(Matrix.toString(Matrix.balancedEigenValue(a)));
System.out.println("Polinomun kokleri : (poly_roots) ");
System.out.println(Matrix.toString(Matrix.poly_roots(p)));
System.out.println("Polinomun kokleri : (poly_rootsC) ");
System.out.println(Matrix.toString(Matrix.poly_rootsC(p)));
}
}
cözüm :
Orijinal Polinom :
-1.0 8.0 -19.0 12.0
Karekteristik matris :
8.0 -19.0 12.0
1.0 0.0 0.0
0.0 1.0 0.0
Eigen Value (sinir deger) :
3.9999999999999964 3.0000000000000067 0.9999999999999991
0.0 0.0 0.0
Eigen Value (sinir deger) balancedeigenValue metodu :
4.000000000000007 2.999999999999993 1.0000000000000004
0.0 0.0 0.0
Polinomun kokleri : (poly_roots)
4.000000000000007 2.999999999999993 1.0000000000000004
0.0 0.0 0.0
Polinomun kokleri : (poly_rootsC)
4.000000000000007 2.999999999999993 1.0000000000000004
Determinant ,
Problem 12.5
Matrix5.java
, matris çarpım, inverse ve determinant işlemleri
import java.io.*;
import Numeric;
class Matrix5
{
public static void main(String args[]) throws IOException
{
//calculation of an inverse matrix
double a[][];
a=new double[2][2];
a[0][0]=1;
a[0][1]=2;
a[1][0]=3;
a[1][1]=4;
System.out.println("Matrix : ");
System.out.println(Matrix.toString(a));
System.out.println("Inverse Matrix : ");
System.out.println(Matrix.toString(Matrix.inv(a)));
System.out.println("Matrix * Inverse Matrix :");
System.out.println(Matrix.toString(Matrix.multiply(a,Matrix.inv(a))));
System.out.println("Determinant = "+Matrix.det(a));
}
}
cözümü :
Matrix :
1.0 2.0
3.0 4.0
Inverse Matrix :
-2.0 1.0
1.5 -0.5
Matrix * Inverse Matrix :
1.0 0.0
0.0 1.0
Determinant = -2.0
Eigenvalue, EigenVektör : Burada verilen matris simetrik bir matristir. sınır değer vektörü QL metodu kullanılarak hesaplanmıştır. Simetrik matris olduğundan sınır değerler gerçek değerler olarak elde edilmiştir.
Problem 12.5 Matrix6.java
, eigenvalue (sınır-değer) problemi
import java.io.*;
import Numeric;
class Matrix6
{
public static void main(String args[]) throws IOException
{
//Eigen Value calculations using QL method
double a[][]={{3,-1,0},{-1,2,-1},{-0,-1,3}};
System.out.println("Original Matrix : ");
System.out.println(Matrix.toString(a));
System.out.println("eigen values and Eigen Vectors : ");
System.out.println(Matrix.toString(Matrix.eigenQL(a)));
}
}
çözümü
:
Original Matrix
:
3.0 -1.0 0.0
-1.0 2.0 -1.0
0.0 -1.0 3.0
eigen values and Eigen
Vectors :
4.0 2.9999999999999996
1.0000000000000004
1.0 1.0 1.0
-1.0000000000000004 0.0 2.0000000000000004
1.0000000000000007-0.9999999999999997 1.0000000000000004
12.3 NUMERIC SINIFI (SAYISAL ANALIZ PAKETİ)
Numeric sınıfı genel bir sayısal analiz paketi olarak düşünülmüştür. Bu pakette benim kişisel ihtiyaçlarım için yazdığım yada başka dillerde daha önceden yazılmış olan kodlardan adepte ettiğim bazı programlar verilmiştir. Bir paketin genel olarak herhangi bir fonksiyonla kullanılabilmesi için herhangi bir fonksiyonla birlikte kullanılması gerekir. Bu yüzden programlarda fonksiyonlar tanımlanmış abstract sınıflar içinde yer almıştır. Gerçek fonksiyonlar yazılırken abstract sınıflardan türetilmeleri gerekmektedir. Numeric paketindeki ve abstract sınıflar ve metodlar şu şekilde tanımlanmıştır :
Numerik sınıfı ve fonksiyon girdi absract sınıfları ve metotları :
abstract class f_x
tek değişkenli tek fonksiyonları aktarmak için kullanılan abstract sınıftır.
Bu sınıfla kullanılabilecek fonksiyona örnek olarak f=x*x fonksiyonu
verilebilir.
Buradaki abstract fonksiyonun tanımı :
abstract
double func(double x);
şeklindedir. Örnek verecek olursak :
class f1 extends f_x
{ double func(double x) { return x*x; } }
Bu sınıf dinamik olarak tanımlandığından kullanılacağında tanımlanması gerekir.
double y=Numeric.derivative(new f1(),x);
veya
f1 xkare=new f1();
double y=Numeric.derivative(xkare,x);
şeklinde çağırılır.
abstract class f_xi
birden fazla fonksiyon birden fazla değişken ile kullanılacaksa f_xi
abstract sınıfı kullanılır. Bu sınıfta tanımlanan func metodu denklem_referansı
indisinin değerine göre tek bir fonksiyon girdisi verir.
örnek
f[0] = x[0]+sin(x[1]);
f[1] = x[0]*x[0]-x[1] ;
func(x,1) f[1]değerini
verir
func(x,0) f[0]değerini verir
abstract metodun tanımı
:
abstract double func(double
x[],int equation_ref);
şeklindedir.
Örnek problemin java kodu eşdeğeri :
class f2 extends f_xi
{ double func(double x[],int x_ref) {
double a=0;
switch(x_ref) { {case 0: a = x[0]+sin(x[1]); break;}
{case 1: a = x[0]*x[0]-x[1]; break;}}
return a; } }
şeklinde verilebilir. Bu funksiyonun asıl programda çağırılması :
double x[]={1,2};
double y=Numeric.derivative(new f2(),x,0);
veya
double x[]={1,2};
f2 xkare=new f2();
double y=Numeric.derivative(xkare,x,0);
şeklinde gerçekleştirilir.
abstract class fi_xi
birden fazla fonksiyon birden fazla değişken ile kullanılacaksa fi_xi
abstract sınıfı kullanılır. Bu sınıfta tanımlanan func metodu tüm fonksiyon
değerlerini vektör olarak iletir.
örnek
f[0]=x[0]+sin(x[1]);
f[1]=x[0]*x[0]-x[1];
func(x) f[0] ve f[1]
değerlerini verir.
abstract double[]
func(double x[]);
Örnek problemin java kodu eşdeğeri :
class f3 extends fi_xi
{ double[] func(double x[])
{ double a[]=new double[2];
a[0] = x[0]+sin(x[1]);
a[1] = x[0]*x[0]-x[1];
return a; }
}
şeklinde verilebilir. Bu funksiyonun asıl programda çağırılması :
double x[]={1,2};
double y[]=Numeric.derivative(new f3(),x);
veya
double x[]={1,2};
f3 fx=new f3();
double y[]=Numeric.derivative(fx,x);
şeklinde gerçekleştirilir.
abstract class
fij_xi
{
iki boyutlu fonksiyon matrisi ve bir boyutlu değişken vectörü verilmişse
fij_xi sınıfını kullanabiliriz. Bu tür fonksiyonlar özellikle iki boyutlu
fonksiyon setinin türev fonksiyonlarını hesaplarken oluşabilir. Çıktı
matris olarak tüm fonksiyonları verir.
örnek :
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]
f[0][0], f[0][1]
f[1][0], f[1][1]
değerlerin iki boyutlu
matrix olarak verir.
abstract double[][] func(double x[]);
örneğin java kodu eşdeğeri :
class f4 extends fij_xi
{
double[][] func(double x[])
{
double b[][];
f=new double[2][2];
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];
return f;
}
}
metod ve sınıfın java programında kullanılışı :
double x[]={1,2};
double y[][]=Numeric.derivative(new f4(),x);
veya
double x[]={1,2};
f4 fx=new f4();
double y[][]=Numeric.derivative(fx,x);
public class Numeric
Numeric sınıfındaki metodlar genelde statik metodlardır. Bu yüzden bağımsız olarak kullanılabilirler. Şimdi bu metodları daha detaylı inceleyelim :
Türev
: bu metotlar gurubu yukarıda
tanımlanan double sayı, vektor ve matris olarak verilen bir fonksiyonun
türevlerini alır.
public static double
derivative(f_x f_deriv,double
x)
df/dx
tek fonksiyon,
tek değişken türevi, örnek fonksiyon f=x*x
fonksiyon abstract sınıfı : f_x
türevin alınacağı nokta : double x
public static double
derivative(f_xi f_deriv,int equation_ref,
double x[],int x_ref)
dfi/ dxj
birden fazla fonksiyon birden fazla değişken türevi
denklem equation_ref indisinin ve x_ref x indisinindeğerine göre
sadece tanımlanan fonksiyonun tanımlanan
x_ref indisine göre olan türevini verir
örnek fonksiyon
f[0] = x[0]+sin(x[1]);
f[1] = x[0]*x[0]-x[1] ;
fonksiyon abstract sınıfı : f_xi
yukarıdaki örnek
için türevin alınacağı nokta : double x[] ile belirlenmiştir
public static double
derivative(fi_xi f_deriv,
int equation_ref,double x[],int x_ref)
public static double
derivative(fi_xi f_deriv,
int equation_ref,double x[],int x_ref)
dfi/ dxj
birden fazla fonksiyon birden fazla değişken türevi
denklem equation_ref indisinin ve x_ref x indisinindeğerine göre
sadece tanımlanan fonksiyonun tanımlanan
x_ref indisine göre olan türevini verir
örnek fonksiyon
f[0] = x[0]+sin(x[1]);
f[1] = x[0]*x[0]-x[1] ;
fonksiyon abstract sınıfı : f_xi
yukarıdaki örnek
için türevin alınacağı nokta : double x[] ile belirlenmiştir
not1 : bir üsteki türev formülüyle bu türev formülünün temel farkı, burada tanımlanan fonksiyonun tüm vektörü çıktı olarak vermesidir. sonuçlar her ikisinde de aynıdır.
not2: denklemlerin hassasiyeti metotların içerisinde verilen h0 ve n değişkenlerinin değiştirilmesiyle arttırılabilir. Şu andaki değerleri optimize edilmiştir ve denenen tüm fonksiyonlar için hassas değerler vermiştir. Nümerik türevlerle çözülmesi oldukça zor olan boyutlu newton-raphson tekniğinde bile başarılı sonuçlar vermiştir. Daha küçük değerler kullanırken önce türevini bildiğiniz fonksiyonlarla sonuçları kontrol etmeniz tavsiye edilir.
İntegral:
Sayısal olarak integral almak için çok çeşitli integrasyon metotları
bulunmaktadır. Burada verilen metodlarda ana kriteria sonuçların hassaslığıdır.
Denklemler nümerik olarak doğru sonuç çıkarması oldukça zor olan siyah cisim
radyasyon fonksiyonu ve istatistik fonksiyonlarda bile tam hassaslık vermiştir.
integral 30 noktalı Gauss - Legendre integral formülü kullanmaktadır. trap,
trapezoidal integral formülü romberg integrasyonunun temelini teşkil etmek
için hazırlanmıştır, tekbaşına da kullanılabilir. Çok hassasiyet gereken
uygulamalarınız için kullanmanızı tavsiye etmiyoruz. Romberg integrasyonu
iterativ bir integrasyon tekniğidir. Trapezoidal integrasyonla hesaplanan
kök değerlerini alarak iterasyonla hata miktarını azaltır.
public static double integral (f_x f_xnt,double a,double b)
f(x) fonksiyonunun double a ve double b sınırları arasında integrali, Gauss-Legendre metodu
örnek fonksiyon
f=x*x
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.integral(new fx(),a,b);
public static void trap
(f_x ff,double a,double h, int j,int m,double R[][])
f(x) fonksiyonunun a ve b sınırları arasında integrali, trapezoidal metodu
örnek fonksiyon
f=x*x
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.trap(new fx(),a,h,j,m,R);
public static double integral_romberg
(f_x ff,double a,double b)
f(x) fonksiyonunun a ve b sınırları arasında integrali, Romberg integrasyon metodu
örnek fonksiyon
f=x*x
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.integral_romberg(new fx(),a,b);
Non-linear denklem
kökleri : non linear denklem
veya denklem sistemlerinin çözümü nümerik analizde en çok karşımıza çıkan
problemlerden biridir.
Genelde tüm metodlar iteratif yöntemler kullanırlar. Bu metodlar içinde
en verimli olarak hesap yapan metodun newton metodu olduğu söylenebilir,
fakat newton metodu fonksiyonun türevinin de hesaplanmasını gerekli kılar.
buradaki metodların bir kısmında türev girdi olarak alınırken, bir kısmında
da yine nümerik metodlarla hesaplanmıştır. denklemin ikinci türevini de
göz önüne alan newton_2nd_derivative gibi metodlar da verilmiştir.
ayrıca secant, bisection gibi türev hesabı gerektirmeyen metodlar da mevcuttur.
Bisection metodu temelde bir arama metodudur ve verimi en düşük metoddur,
fakat verilen bölgeyi tam olarak taradığı ve bir kök varsa muhakkak bulabildiği
için verimsiz bir metod için oldukça sık tercih gören bir yöntemdir.
public static double newton (double x,f_x y,f_x dy)
y(x) fonksiyonunun kökleri, dy(x)/dx fonksiyonu da verilmektedir.
örnek fonksiyon y=x*x-2
dy=2*x
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.newton(x,new fx(),new dfx());
public static double newton(double x,f_x y)
y(x) fonksiyonunun kökleri, dy(x)/dx fonksiyonu nümerik olarak
hesaplanmaktadır
örnek fonksiyon y=x*x-2
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.newton(x,new fx());
public static double newton_2nd_derivative(double x,f_x y,f_x dy)
y(x) fonksiyonunun kökleri, 2inci türev formülü, dy(x)/dx fonksiyonu
da verilmektedir.
örnek fonksiyon y=x*x-2
dy=2*x
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.newton_2nd _derivative (x,new fx(),new dfx(),new dfx() );
public static double
newton_2nd_derivative(double x,f_x y)
y(x) fonksiyonunun kökleri, 2nci türev formülü, dy(x)/dx fonksiyonu
nümerik olarak hesaplanmaktadır
örnek fonksiyon y=x*x-2
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.newton_2nd _derivative (x,new fx(),new dfx());
public static double secant_2nd_derivative(double x,f_x y)
y(x) fonksiyonunun kökleri, 2nci türev formülü, secand metodu kullanılmıştır (türev formülü gerekmez)
örnek fonksiyon y=x*x-2
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric.secant_2nd _derivative (x,new fx(),new dfx());
public static double bisection
(double xl,double xu,f_x y)
y(x) fonksiyonunun kökleri, 2nci türev formülü, bisection metodu kullanılmıştır
(türev formülü gerekmez, temel olarak bir arama metodudur.)
örnek fonksiyon y=x*x-2
fonksiyon abstract sınıfı : f_x
örnek çağırma : double y=Numeric. bisection (xl,xu,new fx());
public static double[] newton( double x[],fi_xi y,fij_xi dy)
yi(xi) fonksiyon sisteminin kökleri kökleri,
dyi(x)/dxj fonksiyonu da matrix fonksiyonu olarak
verilmektedir.
örnek fonksiyon
y1=2*x1*x1*x2-3*x1
y2=2*x1+2*x2
örnek türev fonksiyonu
dy1/dx1=4*x1*x2-3 dy1/dx2=2*x1
dy2/dx1=2 dy2/dx2=2
fonksiyon abstract sınıfı : fi_xi
türev fonksiyonu abstract sınıfı . fi_xij
örnek çağırma : double y[]=Numeric. newton(x,new fixi(),new dfixij());
public static double[] newton( double x[],fij_xi y)
yi (xi) fonksiyon sisteminin kökleri kökleri, dyi (x)/dx j fonksiyonu nümerik olarak hesaplanmaktadır, girdi değildir.
örnek fonksiyon
y1=2*x1*x1*x2-3*x1
y2=2*x1+2*x2
fonksiyon abstract sınıfı : fi_xi
örnek çağırma : double y[]=Numeric. newton(x,new fixi());
Least Square (en
küçük kareler metodu) polinom uydurma
public static double[]
poly_least_square (double xi[],double yi[],int n)
polynom least square (en küçük kareler metodu) polinom uydurma,
çıktı olarak polinomun katsayılarını verir
örnek çağırma : double y[]=Numeric. poly_least_square (xi,yi,n);
public static double f_least_square(double e[],double x)
polynom least square (en küçük kareler metodu) polinom uydurma,
polinomun katsayılarını kullanarak fonksiyon değerini hesaplar
örnek çağırma : double y=Numeric.
f_least_square
(e,x);
public static double error(double x[],double y[],double e[])
polynom least square (en küçük kareler metodu) polinom uydurma,
polinomun katsayılarını kullanarak fonksiyon değerini hesaplar ve toplam
hata miktarını bulur.
örnek çağırma
: double y=Numeric.
error(x,y,e);
public static double[][]
orthogonal_polynomial_least_square (double xi[],double fi[],int m)
ortagonal polynom least square (en küçük kareler metodu) polinom
uydurma, çıktı olarak polinomun
katsayılarını verir. hesaplama yöntemi matrislere dayanmadığından matris çözümlerine dayanan hataları
yapmaz. yalnız polinom formülü daha kompleksdir
örnek çağırma : double y[]=Numeric. Orthogonal_polynomial_least_square (xi,yi,n);
public static double f_orthogonal_least_square(double e[],double x)
ortogonal polynom least square (en küçük kareler metodu) polinom uydurma, polinomun katsayılarını
kullanarak fonksiyon değerini hesaplar
örnek
çağırma : double y=Numeric.
f_orthogonal_least_square
(e,x);
Diferansiyel denklemler
public static double[] RK4(f_xi fp,double x0,double xn,double f0,int N)
4üncü dereceden Runge-Kutta diferansiyal denklem çözüm metodu,
tek bir denklemi çözmek için
örnek
çağırma : double y[]=RK4 (new fxi(),x0,xn,N);
public static double[][] RK4(f_xi fp,double x0,double xn,double f0[],int N)
4üncü dereceden Runge-Kutta diferansiyal denklem sistemi çözüm metodu, bir denklemin sistemini veya
daha yüksek
dereceden diferansiyel denklemleri çözmek için kulanılır.
örnek çağırma : double y[][]=RK4 (new fxi(),x0,xn,f0,N);
public static double[][] RKF45(f_xi fp,double x0,double xn,double f0,int N) throws IOException
4-5 incü dereceden Runge-Kutta-fehlenberg diferansiyal denklem sistemi çözüm metodu, bir denklemin
sistemini
veya daha yüksek dereceden diferansiyel denklemleri çözmek için kulanılır.
örnek çağırma : double y[][]=RKF45 (new fxi(),x0,xn,f0,N);
Problem 12.6 Numeric
sınıfı listesi :
Numeric.java
//=============================================================
// Numerical Analysis Package in Java
// Numerical analysis calculation class Numeric
// Author : Dr. Turhan Coban
//=============================================================
import java.io.*;
import Matrix;
abstract class f_x
{
//single function single independent variable
// example f=x*x
abstract double func(double x);
}
abstract class f_xi
{
// multifunction multi independent variable
// a single value is returned indiced to equation_ref
// example f[0]=x[0]+sin(x[1])
// f[1]=x[0]*x[0]-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[],int equation_ref);
}
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[]);
}
public class Numeric
{
//This is a library of numerical methods
//
public static double derivative(f_x f_deriv,double x)
{
// This method calculates derivatives of a simple function
// accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_x
double h0=0.0256;
// accuracy can be change by adjusting h0 and n
int i,m;
int n=7;
//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++)
{
T[i][0]=( f_deriv.func(x + h[i]) - f_deriv.func( x - h[i]))/(2.0*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
derivative(f_xi f_deriv,int equation_ref,
double x[],int x_ref)
{
// This method calculates derivative of a function selected from a set
of
// functions. Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
// df(equation_ref)/dx(x_ref)
double h0=0.0256;
int i,m;
int n=7;
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];
T[i][0]=( f_deriv.func(x1,equation_ref) -
f_deriv.func(x2,equation_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 static double
derivative(fi_xi f_deriv,
int equation_ref,double x[],int x_ref)
{
// This method calculates derivative of a function selected from a set
of
// functions. Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
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_deriv.func(x1);
f2=f_deriv.func(x2);
T[i][0]=( f1[equation_ref] - f2[equation_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 static double
integral(f_x f_xnt,double a,double b)
{
//integral of a function by using gauss-legendre quadrature
//
double s[],w[];
int i;
s=new double[30];
w=new double[30];
s[ 0] = .15532579626752470000E-02;
s[ 1] = .81659383601264120000E-02;
s[ 2] = .19989067515846230000E-01;
s[ 3] = .36899976285362850000E-01;
s[ 4] = .58719732103973630000E-01;
s[ 5] = .85217118808615820000E-01;
s[ 6] = .11611128394758690000E+00;
s[ 7] = .15107475260334210000E+00;
s[ 8] = .18973690850537860000E+00;
s[ 9] = .23168792592899010000E+00;
s[10] = .27648311523095540000E+00;
s[11] = .32364763723456090000E+00;
s[12] = .37268153691605510000E+00;
s[13] = .42306504319570830000E+00;
s[14] = .47426407872234120000E+00;
s[15] = .52573592127765890000E+00;
s[16] = .57693495680429170000E+00;
s[17] = .62731846308394490000E+00;
s[18] = .67635236276543910000E+00;
s[19] = .72351688476904450000E+00;
s[20] = .76831207407100990000E+00;
s[21] = .81026309149462140000E+00;
s[22] = .84892524739665800000E+00;
s[23] = .88388871605241310000E+00;
s[24] = .91478288119138420000E+00;
s[25] = .94128026789602640000E+00;
s[26] = .96310002371463720000E+00;
s[27] = .98001093248415370000E+00;
s[28] = .99183406163987350000E+00;
s[29] = .99844674203732480000E+00;
w[ 0] = .39840962480827790000E-02;
w[ 1] = .92332341555455000000E-02;
w[ 2] = .14392353941661670000E-01;
w[ 3] = .19399596284813530000E-01;
w[ 4] = .24201336415292590000E-01;
w[ 5] = .28746578108808720000E-01;
w[ 6] = .32987114941090080000E-01;
w[ 7] = .36877987368852570000E-01;
w[ 8] = .40377947614710090000E-01;
w[ 9] = .43449893600541500000E-01;
w[10] = .46061261118893050000E-01;
w[11] = .48184368587322120000E-01;
w[12] = .49796710293397640000E-01;
w[13] = .50881194874202750000E-01;
w[14] = .51426326446779420000E-01;
w[15] = .51426326446779420000E-01;
w[16] = .50881194874202750000E-01;
w[17] = .49796710293397640000E-01;
w[18] = .48184368587322120000E-01;
w[19] = .46061261118893050000E-01;
w[20] = .43449893600541500000E-01;
w[21] = .40377947614710090000E-01;
w[22] = .36877987368852570000E-01;
w[23] = .32987114941090080000E-01;
w[24] = .28746578108808720000E-01;
w[25] = .24201336415292590000E-01;
w[26] = .19399596284813530000E-01;
w[27] = .14392353941661670000E-01;
w[28] = .92332341555455000000E-02;
w[29] = .39840962480827790000E-02;
int n=30;
double z=0;
double x,y;
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*s[i];
y=f_xnt.func(x);
z+=(b-a)/2*w[i]*y;
}
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*(-s[i]);
y=f_xnt.func(x);
z+=(b-a)/2.0*w[i]*y;
}
return z;
}
public static void
trap(f_x ff,double a,double h, int j,int m,double R[][])
{
// successive trapezoidal integration rule
// this program will be utilised in romberg integration
double sum=0;
int p;
for(p=1;p<=m;p++)
{ sum+=ff.func(a+h*(2*p-1)); }
R[j][0]=R[j-1][0]/2.0+h*sum;
}
public static double
integral_romberg(f_x ff,double a,double b)
{
//romberg integration
int n=8;//increase n to increase accuracy
double R[][];
R=new double[n+1][n+1];
int m=1;
double h=b-a;
double close=1;
double tol=1e-40;
int j=0,k=0;
double ret=0;
R[0][0]=h/2.0*(ff.func(a)+ff.func(b));
do
{
j++;
h/=2.0;
Numeric.trap(ff,a,h,j,m,R);
m*=2;
for(k=1;k<=j;k++)
{
R[j][k]=R[j][k-1]+(R[j][k-1]-R[j-1][k-1])/(Math.pow(4,k)-1.0);
close=Math.abs(R[j-1][j-1]-R[j][j]);
ret=R[j][k];
}
}while(closetol && j<n);
return ret;
}
//==================================================
//Finding Non-linear Roots of Functions
public static double
newton(double x,f_x y,f_x dy)
{
// Newton - Raphson
method for single equation
// with single variable
// required function and its derivative
int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=dy.func(x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
public static double
newton(double x,f_x y)
{
// Newton - Raphson method for single equation
// required function only derivative is calculated
// numerically by using method derivative
int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=Numeric.derivative(y,x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
public static double
newton_2nd_derivative(double
x,f_x y,f_x dy)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
//function and first order derivative is required
int nmax=500;
double dx=1e-15;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=dy.func(x);
dfxm=dy.func(x-dx);
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
public static double
newton_2nd_derivative(double x,f_x y)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function is required, first order derivative calculated
// numerically by derivative method.
int nmax=500;
double dx=1e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=Numeric.derivative(y,x);
dfxm=Numeric.derivative(y,(x-dx));
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
public static double
secant_2nd_derivative(double
x,f_x y)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function should be supplied
int nmax=500;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
fxm1=y.func(x-dx1);
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
public static double
bisection(double xl,double xu,f_x y)
{
//bisection method to find roots of a function y.func(x)
//function should be supplied
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter : number of iteration
double test;
double xr=0;
double es,ea;
double fxl,fxr;
int maxit=500,iter=0;
es=0.000001;
ea=1.1*es;
while((eaes)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
{ ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl= y.func(xl);
fxr= y.func(xr);
test= fxl*fxr;
if(test==0.0) ea=0;
else if(test<0.0) xu=xr;
else
{
xl=xr;
}
} //end of while
return xr;
}
public static double[]
newton( double x[],fi_xi y,fij_xi dy)
{
//solution of nonlinear system of equations
//by using newton raphson method.
//ti :weigth factor
//x independent value vector:
//y dependent function vector
//dy derivative of dependent function matrix
double ti=1.0;
int i;
int nmax=500;
double tolerance=1.0e-30;
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(y.func(x),dy.func(x)),-ti);
x=Matrix.add(x,b);
}
if(i = nmax) System.out.println("warning maximum number"+
" of iterations results may not be valid");
return x;
}
public static double[]
newton( double x[],fi_xi y)
{
//solution of nonlinear system of equations
//by using newton raphson method.
//this function does not require derivatives
//it is utilised method derivative to calculate derivatives
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 )
{
for(ii=0;ii<n;ii++)
{
for(jj=0;jj<n;jj++)
{
dy[ii][jj]=Numeric.derivative(y,ii,x,jj);
}
}
b=Matrix.multiply(Matrix.divide(y.func(x),dy),-ti);
x=Matrix.add(x,b);
}
if(i = nmax) System.out.println("warning maximum number of iterations"+
" results may not be valid"
);
return x;
}
//========= least square curve fitting methods ==========
//-------- polynomial least square curve fitting --------
public static double[]
poly_least_square
(double xi[],double yi[],int n)
//polynomial least square curve fitting
//variables xi,yi vector of data points
//degree of curve fitting
{
//xi vector of independent variable
//yi vector of dependent variable
//n : degree of curve fitting
int l=xi.length;
int i,j,k;
int np1=n+1;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double X[];
X=new double[np1];
for(i=0;i<n+1;i++)
{
for(j=0;j<n+1;j++)
{
if(i==0 && j==0)
A[i][j]=l;
else
for(k=0;k<l;k++)
A[i][j] += Math.pow(xi[k],(i+j));
}
for(k=0;k<l;k++)
{
if(i==0)
B[i]+= yi[k];
else
B[i] += Math.pow(xi[k],i)*yi[k];
}
}
X=Matrix.divide(B,A);
//X=B/A;
double max=0;
for(i=0;i<n+1;i++)
if(Math.abs(X[i]) max) max = Math.abs(X[i]);
for(i=0;i<n+1;i++)
if((Math.abs(X[i]/max) 0) && (Math.abs(X[i]/max) < 1.0e-11))
X[i]=0;
return X;
}
public static double
f_least_square(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=e.length;
double ff;
if(n!=0.0)
{ ff=e[n-1];
for(int i=n-2;i=0;i--)
{ ff=ff*x+e[i]; }
}
else
ff=0;
return ff;
}
public static double
error(double x[],double y[],double e[])
{
//calculates absolute square root error of a least square approach
double n=x.length;
int k;
double total=0;
for(k=0;k<n;k++)
{
total+=(y[k]-f_least_square(e,x[k]))*(y[k]-f_least_square(e,x[k]));
}
total=Math.sqrt(total);
return total;
}
//-------End least square
methods --------------------
//-------Orthogonal polynomial least square curve fitting ---
public static double[][]
orthogonal_polynomial_least_square
(double xi[],double fi[],int m)
//orthogonal polynomial least square curve fitting
//this method do not require any matrix solution
//variables xi,fi vector of data points wi weight functions
//m degree of curve fitting
{
int i,j,k;
int n=xi.length;
int mp2=n+2;
int mp1=n+1;
//matrix<double p(m+2,n,0);
double p[][];
p=new double[mp2][n];
//vector<double gamma(m+1,0);
double gamma[];
gamma=new double[mp1];
//vector<double beta(m+1,0);
double beta[];
beta=new double[mp1];
//vector<double omega(m+1,0);
double omega[];
omega=new double[mp1];
//vector<double alpha(m+1,0);
double alpha[];
alpha=new double[mp1];
//vector<double b(m+1,0);
double b[];
b=new double[mp1];
//vector<double wi(n,1.);
double wi[];
wi=new double[n];
//matrix<double a(3,m+1);
double a[][];
a=new double[3][mp1];
for(i=0;i<n;i++)
{
p[1][i]=1.0;
p[0][i]=0.0;
}
gamma[0]=0;
for(i=0;i<n;i++)
{
gamma[0]+=wi[i];
}
beta[0]=0.0;
for(j=0;j<m+1;j++)
{
omega[j]=0;
for(i=0;i<n;i++)
{
omega[j]+=wi[i]*fi[i]*p[j+1][i];
}
b[j]=omega[j]/gamma[j];
if( j != m)
{
alpha[j+1]=0;
for(i=0;i<n;i++)
{
alpha[j+1]+=wi[i]*xi[i]*p[j+1][i]*p[j+1][i]/gamma[j];
}
for(i=0;i<n;i++)
{
p[j+2][i]=(xi[i]-alpha[j+1])*p[j+1][i]-beta[j]*p[j][i];
}
gamma[j+1]=0;
for(i=0;i<n;i++)
{
gamma[j+1]+=wi[i]*p[j+2][i]*p[j+2][i];
}
beta[j+1]=gamma[j+1]/gamma[j];
}
}//end of j
for(j=0;j<m+1;j++)
{
a[0][j]=b[j];
a[1][j]=alpha[j];
a[2][j]=beta[j];
}
return a;
}
public static double
f_orhogonal_polynomial(double a[][],double x)
{
double yy=0;
int k;
int m=a[0].length-1;
int mp2=m+2;
double q[];
q=new double[mp2];
//vector<double q(m+2,0.0);
for(k=m-1;k=0;k--)
{
q[k]=a[0][k]+(x-a[1][k+1])*q[k+1]-a[2][k+1]*q[k+2];
yy=q[k];
}
return yy;
}
//minimization of a
function
public static double amotry(double p[][], double y[], double psum[],
int ndim,
f_xi ff, int ihi, double fac)
{
int j;
double fac1,fac2,ytry;
double ptry[]=new double[ndim];
fac1=(1.0-fac)/ndim;
fac2=fac1-fac;
for (j=1;j<=ndim;j++) ptry[j-1]=psum[j-1]*fac1-p[ihi-1][j-1]*fac2;
ytry=ff.func(ptry,0);
if (ytry < y[ihi-1]) {
y[ihi-1]=ytry;
for (j=1;j<=ndim;j++) {
psum[j-1] += ptry[j-1]-p[ihi-1][j-1];
p[ihi-1][j-1]=ptry[j-1];
}
}
return ytry;
}
public static void
amoeba(double p[][], double y[], int ndim, double ftol,
f_xi ff, int nfunk)
{
int i,ihi,ilo,inhi,j,mpts=ndim+1;
int NMAX=5000;
double rtol,sum,swap,ysave,ytry;
double psum[]=new double[ndim];
nfunk=0;
for (j=1;j<=ndim;j++)
{
for (sum=0.0,i=1;i<=mpts;i++) sum += p[i-1][j-1];
psum[j-1]=sum;
}
for (;;) {
ilo=1;
//ihi = y[1]y[2] ? (inhi=2,1) : (inhi=1,2);
if(y[0]y[1])
{
inhi=2;
ihi=1;
}
else
{
inhi=1;
ihi=2;
}
for (i=1;i<=mpts;i++) {
if (y[i-1] <= y[ilo-1]) ilo=i;
if (y[i-1] y[ihi-1]) {
inhi=ihi;
ihi=i;
} else if (y[i-1] y[inhi-1] && i != ihi) inhi=i;
}
rtol=2.0*Math.abs(y[ihi-1]-y[ilo-1])
/(Math.abs(y[ihi-1])+Math.abs(y[ilo-1]));
if (rtol < ftol) {
swap=y[0];
y[0]=y[ilo-1];
y[ilo-1]=swap;
for (i=1;i<=ndim;i++)
{
swap=p[0][i-1];
p[0][i-1]=p[ilo-1][i-1];
p[ilo-1][i-1]=swap;
}
break;
}
if (nfunk = NMAX) System.out.println("NMAX exceeded");
nfunk += 2;
ytry=amotry(p,y,psum,ndim,ff,ihi,-1.0);
if (ytry <= y[ilo-1])
ytry=amotry(p,y,psum,ndim,ff,ihi,2.0);
else if (ytry = y[inhi-1]) {
ysave=y[ihi-1];
ytry=amotry(p,y,psum,ndim,ff,ihi,0.5);
if (ytry = ysave) {
for (i=1;i<=mpts;i++) {
if (i != ilo) {
for (j=1;j<=ndim;j++)
psum[j-1]=0.5*(p[i-1][j-1]+p[ilo-1][j-1]);
p[i-1][j-1]=psum[j-1];
y[i]=ff.func(psum,0);
}
}
nfunk += ndim;
for (j=1;j<=ndim;j++)
{
for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];
psum[j]=sum;
}
}
} else --nfunk;
}
}
// differential equation
solutions
public static double[]
RK4(f_xi fp,double x0,double xn,double f0,int N)
{
//4th order Runge Kutta Method
//fp :given derivative function df/dx(f,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;
double fi[]=new double[N+1];
double xi[]=new double[2];
double k[];
k=new double[4];
int i;
double x;
fi[0]=f0;
for(x=x0,i=0;x<xn;x+=h,i++)
{
xi[0]=x;
xi[1]=fi[i];
k[0]=h*fp.func(xi,0);
xi[0]=x+h/2.0;
xi[1]=fi[i]+k[0]/2;
k[1]=h*fp.func(xi,0);
xi[1]=fi[i]+k[1]/2.0;
k[2]=h*fp.func(xi,0);
xi[0]=x+h;
xi[1]=fi[i]+k[2];
k[3]=h*fp.func(xi,0);
fi[i+1]=fi[i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
}
return fi;
}
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];
for(x=x0,i=0;x<xn;x+=h,i++)
{
for(j=1;j<=M;j++)
{
xi[0]=x;
xi[1]=fi[j-1][i];
k[0]=h*fp.func(xi,j-1);
xi[0]=x+h/2.0;
xi[1]=fi[j-1][i]+k[0]/2;
k[1]=h*fp.func(xi,j-1);
xi[1]=fi[j-1][i]+k[1]/2.0;
k[2]=h*fp.func(xi,j-1);
xi[0]=x+h;
xi[1]=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;
}
}
return fi;
}
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(s1.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;
}
//End of Class Numeric
}
12.4 NUMERIC SINIFI (SAYISAL ANALIZ PAKETİ) ÖRNEK PROBLEMLERİ
Non linear denklem çözümleri :
Problem
12.7 non lineer denklem sistemi çözümü örneği ,
Matrix2.java
import java.io.*;
import Numeric;
import complex;
class f1 extends fi_xi
{
double[] func(double[] x)
{
double b[];
b=new double[2];
b[0]=x[0]*x[0]-x[1]-1;
b[1]=(x[0]-2)*(x[0]-2)+(x[1]-0.5)*(x[1]-0.5)-1.0;
return b;
}
}
class df1 extends fij_xi
{
double[][] func(double x[])
{
double b[][];
b=new double[2][2];
b[0][0]=2.0*x[0];
b[0][1]=-1.0;
b[1][0]=2.0*x[0]-4.0;
b[1][1]=2.0*x[1]-4.0;
return b;
}
}
class f2 extends f_x
{
double func(double x)
{
return 4.2+0.45*x+0.0025*Math.pow(x,2.8)-(60+8*x-0.16*x*x);
}
}
class df2 extends f_x
{
double func(double x)
{
return Math.cos(x)-0.5;
}
}
class Matrix2
{
public static void main(String args[]) throws IOException
{
//Solution of zero of functions of one or multivariable equations
f1 b1;
b1=new f1();
df1 db1;
db1=new df1();
f2 b2;
b2=new f2();
df2 db2;
db2=new df2();
double x[];
x=new double[2];
x[0]=0;
x[1]=0.5;
double xx=13;
System.out.println("Roots of equations : ");
System.out.println("Newtons method y=f(x)=0 f(x) and df(x)/dx
is given");
System.out.println("Newton y,dy x= "+Numeric.newton(xx,b2,db2));
System.out.println("");
System.out.println("Newtons method y=f(x)=0 f(x) is given");
System.out.println("Newton y, x= "+Numeric.newton(xx,b2));
System.out.println("");
System.out.println("Newtons method with 2nd derivative y=f(x)=0
");
System.out.println("f(x) and df(x)/dx is given
");
System.out.println("Newton_2nd_derivative y,dy x= "+
Numeric.newton_2nd_derivative(xx,b2,db2));
System.out.println("");
System.out.println("Newtons method with 2nd derivative y=f(x)=0
");
System.out.println("f(x) is given only ");
System.out.println("Newton_2nd_derivative y, x= "+
Numeric.newton_2nd_derivative(xx,b2));
System.out.println("");
System.out.println("Secant method with 2nd derivative y=f(x)=0
");
System.out.println("f(x) is given only ");
System.out.println("Secant_2nd_derivative y, x= "+
Numeric.secant_2nd_derivative(xx,b2));
System.out.println("");
System.out.println("Newtons method for system of equations");
System.out.println("yi=fi(xi)=0,i=0,n ");
System.out.println("fi(xi) and dfi(xi)/dxj,i=1,n,j=1,n is given");
System.out.println(" "+Matrix.toStringT(Numeric.newton(x,b1,db1)));
System.out.println("");
System.out.println("Newtons method for system of equations");
System.out.println("yi=fi(xi)=0,i=0,n ");
System.out.println("fi(xi),i=1,n is given(derivative matrix is notrequired)");
System.out.println(" "+Matrix.toStringT(Numeric.newton(x,b1)));
}
}
Çözüm seti :
Roots of equations :
Newtons method y=f(x)=0 f(x) and df(x)/dx is given
Newton y,dy x= 1.895494267033981
Newtons method y=f(x)=0 f(x) is given
Newton y, x= 1.895494267033981
Newtons method with 2nd derivative y=f(x)=0
f(x) and df(x)/dx is given
Newton_2nd_derivative y,dy x= 1.895494267033981
Newtons method with 2nd derivative y=f(x)=0
f(x) is given only
Newton_2nd_derivative y, x= 1.895494267033981
Secant method with 2nd derivative y=f(x)=0
f(x) is given only
Secant_2nd_derivative y, x= 1.895494267033981
Newtons method for system of equations
yi=fi(xi)=0,i=0,n
fi(xi) and dfi(xi)/dxj,i=1,n,j=1,n is given
1.0673460858066892
0.13922766688686042
Newtons method for system of equations
yi=fi(xi)=0,i=0,n
fi(xi),i=1,n is given(derivative matrix is notrequired)
1.0673460858066897
0.1392276668868614
İntegral -türev : İntegral ve türev örneği
Problem 12.8 integral ve türev örneği, integ1.java
//======================================================
// Numerical Analysis package in java
// example to show utilisation of integration (integral)
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;
import Numeric;
import complex;
class f1 extends f_x
{
double func(double x)
{
return x*x;
}
}
class f2 extends f_x
{
double func(double x)
{
return x;
}
}
class fm1 extends f_xi
{
//multivariable function
double func(double x[],int x_ref)
{
double a=0;
switch(x_ref)
{
case 0: a=(x[0]+x[1]*x[1]+x[2]*2);break;
case 1: a=(x[0]*3+x[1]*x[1]*x[1]+x[2]);break;
case 2: a=(x[0]-5*x[1]*x[1]+x[2]*x[2]*x[2]);break;
}
return a;
}
}
class fm2 extends fi_xi
{
//multivariable function
double[] func(double x[])
{
double a[];
a=new double[3];
a[0]=(x[0]+x[1]*x[1]+x[2]*2);
a[1]=(x[0]*3+x[1]*x[1]*x[1]+x[2]);
a[2]=(x[0]-5*x[1]*x[1]+x[2]*x[2]*x[2]);
return a;
}
}
class integ1
{
public static void main(String args[]) throws IOException
{
f1 b1;
b1=new f1();
f2 b2;
b2=new f2();
fm1 b3;
b3=new fm1();
double x[];
fm2 b4;
b4=new fm2();
x=new double[3];
x[0]=1;
x[1]=2;
x[2]=3;
System.out.println("integral of class f1 : "+
Numeric.integral(b1,0.0,1.0));
System.out.println("romberg integral of class f1 : "+
Numeric.integral_romberg(b1,0.0,1.0));
System.out.println("integral of class f2 : "+
Numeric.integral(b2,0.0,1.0));
System.out.println("Derivative of class f1 : "+
Numeric.derivative(b1,1.0));
System.out.println("Derivative of class f2 : "+
Numeric.derivative(b2,1.0));
System.out.println("Derivative of fm1,0 : "+
Numeric.derivative(b3,1,x,0));
System.out.println("Derivative of fm2,0 : "+
Numeric.derivative(b4,1,x,0));
System.out.println("Derivative of fm1,0 : "+
Numeric.derivative(b3,1,x,0));
System.out.println("Derivative of fm2,0 : "+
Numeric.derivative(b4,1,x,0));
System.in.read();
}
}
çözüm :
integral of class f1
: 0.333333333333329
romberg integral of class f1 : 0.3333333333333333
integral of class f2 : 0.49999999999999384
Derivative of class f1 : 1.9999999999994864
Derivative of class f2 : 0.9999999999998606
Derivative of fm1,0 : 3.0000000000028932
Derivative of fm2,0 : 3.0000000000028932
Derivative of fm1,0 : 3.0000000000028932
Derivative of fm2,0 : 3.0000000000028932
Least Square polinomu
uydurma
Problem 12.9 least square polinom uydurma, least.java
//===========================================
// Polynomial least square curve fitting
// Dr. Turhan Coban
//===========================================
import java.io.*;
import Numeric;
import Text;
class least
{
public static void main(String args[]) throws IOException
{
//polynomial least square curve fitting
double e[]; //coefficients of polynomial least square
e=new double[20];
double xi[];
xi=new double[500];
double yi[];
yi=new double[500];
DataInputStream cin=new DataInputStream(System.in);
System.out.println("name of input file : ");
String in_name=Text.readString(cin);
System.out.println(" ");
DataInputStream fin=new DataInputStream(new FileInputStream(in_name));
PrintStream fout=new PrintStream(new FileOutputStream("out.dat"));
int i=-1;
try {
while(fin != null)
{
i++;
xi[i]=Text.readDouble(fin);
yi[i]=Text.readDouble(fin);
System.out.println(" "+xi[i]+" "+yi[i]);
}
} catch(EOFException e_eof) {System.out.println("end of file");
}
System.out.println("degree of curve fitting : ");
int n=Text.readInt(cin);
e=Numeric.poly_least_square(xi,yi,n);
System.out.println(Matrix.toString(e));
fout.println(Matrix.toString(e));
System.out.println("least square error : "+Numeric.error(xi,yi,e));
System.out.println("would you like to send curve fitted data to a data
file (y/n) ");
char answer=Text.readChar(cin);
if(answer=='y')
{
System.out.println("name of output file : ");
String out_name=Text.readString(cin);
PrintStream fout1=new PrintStream(new FileOutputStream(out_name));
double xmin,xmax;
System.out.println("minimum x : ");
xmin=Text.readDouble(cin);
System.out.println("maximum x : ");
xmax=Text.readDouble(cin);
System.out.println("number of data points : ");
int n_data=Text.readInt(cin);
double x;
for(i=0;i<n_data;i++)
{
x=xmin+(xmax-xmin)/(n_data-1.0)*i;
System.out.println(" "+x+" "+Numeric.f_least_square(e,x));
fout1.println(" "+x+" "+Numeric.f_least_square(e,x));
} //end of for
}//end of if
}
}
least square polinom uydurma işlevini biraz daha kullanıcıya cazip hale getirmek için grafik formatındaki leastSW.java programı geliştirilmiştir.
Problem 12.10 least square polinom uydurma, leastSW.java jFrame grafik programı
//======================================================
// Numerical Analysis package in java
// En küçük kareler metodu polinom uydurma
// Dr. Turhan Coban
// =====================================================
/*
* Java 2 Swing JFrame version.
*/
import java.lang.Integer;
import java.awt.*;
import java.awt.event.*;
import java.awt.font.*;
import java.awt.geom.*;
import java.awt.image.*;
import java.awt.event.ActionListener;
import java.awt.event.ActionEvent;
import java.awt.event.ItemListener;
import java.awt.event.ItemEvent;
import javax.swing.*;
import javax.swing.border.*;
import Numeric;
import Text;
import java.util.*;
import java.io.*;
import Plot2D;
import PlotShapesSW;
import java.security.*;
public class leastSW extends JFrame implements ItemListener,ActionListener
{
boolean inAnApplet = true;
final static String girdi = "DATA girdi sayfası";
final static String leastcikti = "En kucuk kareler çıktı sayfası";
final static String KONTROLPanel = "PLOT Kontrol sayfası";
final static String PLOTPanel = "Plot sayfası";
Plot2D jta;
StringTokenizer t;
// least değişken seti ********
double ee[];
double xi[];
double yi[];
BufferedReader fin;
File girdiDosyasi;
File ciktiDosyasi;
File plotDosyasi;
int nGirdi;
int nCikti;
double minX,maxX;
int PolinomDerecesi;
// ****** girdi sayfası *******
JLabel promptGirdi;
JTextField inputGirdi;
JLabel promptPD;
JTextField inputPD;
JTextArea outputGirdi;
// ****** çıktı sayfası *******
JLabel promptCikti;
JTextField inputCikti;
JTextArea outputCikti;
JTextArea outputFormul;
JLabel promptP1_1,promptP1_2,promptP1_3,promptP1_4;
JTextField inputP1_1,inputP1_2,inputP1_3,inputP1_4;
// ****** plot kontrol sayfası ****
JLabel promptXmin; // Label prompt in Xmin field
JLabel promptXmax; // Label prompt in Xmax field
JLabel promptYmin; // Label prompt in Ymin field
JLabel promptYmax; // Label prompt in Ymax field
JLabel promptLabel; // Label prompt Plot Label;
JLabel promptXLabel; // Label prompt Plot XLabel;
JLabel promptYLabel; // Label prompt Plot YLabel;
JLabel promptXntic; // Label prompt in Xmin field
JLabel promptYntic; // Label prompt in Xmax field
JLabel promptXgridon; // Label prompt in Ymin field
JLabel promptYgridon; // Label prompt in Ymax field
JTextField inputXmin; // input field Xmin
JTextField inputXmax; // input field Xmax
JTextField inputYmin; // input field Ymin
JTextField inputYmax; // input field Ymax
JTextField inputLabel; // input field Label
JTextField inputXLabel; // input field XLabel
JTextField inputYLabel; // input field YLabel
JTextField inputXntic; // input field xntic
JTextField inputYntic; // input field yntic
JCheckBox inputXgridon; // input field xgridon
JCheckBox inputYgridon; // input field ygridon
JButton printButton;
//===========================================================
public leastSW(String s) throws IOException
{
super(s);
Border b=BorderFactory.createRaisedBevelBorder();
girdiDosyasi=new File("in.txt");
girdiDosyasi=new File(girdiDosyasi.getAbsolutePath());
ciktiDosyasi=new File("out.txt");
ciktiDosyasi=new File(ciktiDosyasi.getAbsolutePath());
plotDosyasi=new File("Plot.txt");
plotDosyasi=new File(plotDosyasi.getAbsolutePath());
oku();
nGirdi=xi.length;
nCikti=nGirdi;
minX=xi[0];
maxX=xi[nGirdi-1];
ee=new double[nGirdi];
PolinomDerecesi=2;
ee=Numeric.poly_least_square(xi,yi,PolinomDerecesi);
Container contentPane = getContentPane();
JTabbedPane tabbedPane = new JTabbedPane();
JPanel pane1 = new JPanel();
JPanel pane2 = new JPanel();
pane2.setLayout(new BorderLayout());
JPanel pane3=new JPanel();
JPanel pane4=new JPanel();
JPanel mpane=new JPanel();
mpane.setLayout(new GridLayout(1,8));
JPanel xpane=new JPanel();
xpane.setLayout(new GridLayout(1,8));
JPanel lpane=new JPanel();
lpane.setLayout(new GridLayout(3,2));
JPanel PanelCikti1=new JPanel();
PanelCikti1.setLayout(new GridLayout(4,2));
JPanel PanelCikti2=new JPanel();
PanelCikti2.setLayout(new BorderLayout());
//JScrollPane.VERTICAL_SCROLLBAR_ALWAYS;
JPanel PanelGirdi1=new JPanel();
PanelGirdi1.setLayout(new GridLayout(2,2));
//***************
promptGirdi=new JLabel("Girdi dosya ismi");
inputGirdi=new JTextField("in.txt");
promptPD=new JLabel("polinom derecesi");
inputPD=new JTextField(""+PolinomDerecesi);
outputGirdi=new JTextArea(yazxiyi(),20,40);
outputGirdi.setLineWrap(true);
outputGirdi.setBorder(b);
promptCikti=new JLabel("Cikti dosya ismi");
inputCikti=new JTextField("out.txt");
promptP1_2=new JLabel("veri sayisi");
promptP1_3=new JLabel("Minimum x : ");
promptP1_4=new JLabel("Maximum x : ");
inputP1_2=new JTextField(""+nCikti);
inputP1_3=new JTextField(""+minX);
inputP1_4=new JTextField(""+maxX);
outputCikti=new JTextArea(yazCikti(),20,40);
outputCikti.setBorder(b);
outputFormul=new JTextArea(yaze(),10,40);
outputFormul.setLineWrap(true);
outputFormul.setAutoscrolls(true);
outputFormul.setBorder(b);
PanelCikti1.add(promptCikti);
PanelCikti1.add(inputCikti);
PanelCikti1.add(promptP1_2);
PanelCikti1.add(inputP1_2);
PanelCikti1.add(promptP1_3);
PanelCikti1.add(inputP1_3);
PanelCikti1.add(promptP1_4);
PanelCikti1.add(inputP1_4);
PanelCikti2.add(PanelCikti1,BorderLayout.NORTH);
PanelCikti2.add(new JScrollPane(outputFormul),BorderLayout.CENTER);
PanelCikti2.add(new JScrollPane(outputCikti),BorderLayout.SOUTH);
PanelGirdi1.add(promptGirdi);
PanelGirdi1.add(inputGirdi);
PanelGirdi1.add(promptPD);
PanelGirdi1.add(inputPD);
pane3.add(PanelGirdi1);
pane3.add(new JScrollPane(outputGirdi));
pane4.add(PanelCikti2);
//***************
promptXmin = new JLabel("Xmin ");
inputXmin = new JTextField(5);
promptXmax = new JLabel("Xmax ");
inputXmax = new JTextField(5);
promptYmin = new JLabel("Ymin ");
inputYmin = new JTextField(5);
promptYmax = new JLabel("Ymax ");
inputYmax = new JTextField(5);
//*******
promptLabel = new JLabel("Plot başlığı : ");
promptXLabel = new JLabel("x ekseni başlığı : ");
promptYLabel = new JLabel("y ekseni başlığı : ");
inputLabel = new JTextField(30);
inputXLabel = new JTextField(30);
inputYLabel = new JTextField(30);
//*******
promptXntic=new JLabel("X tic no");
inputXntic=new JTextField(5);
promptYntic=new JLabel("Y tic no");
inputYntic=new JTextField(5);
promptXgridon=new JLabel("X grid");
inputXgridon=new JCheckBox(" ");
promptYgridon=new JLabel("Y grid");
inputYgridon=new JCheckBox(" ");
//*******
mpane.add(promptXmin);
mpane.add(inputXmin);
mpane.add(promptXmax);
mpane.add(inputXmax);
mpane.add(promptYmin);
mpane.add(inputYmin);
mpane.add(promptYmax);
mpane.add(inputYmax);
pane1.add(mpane,BorderLayout.NORTH);
xpane.add(promptXntic);
xpane.add(inputXntic);
xpane.add(promptYntic);
xpane.add(inputYntic);
xpane.add(promptXgridon);
xpane.add(inputXgridon);
xpane.add(promptYgridon);
xpane.add(inputYgridon);
pane1.add(xpane,BorderLayout.NORTH);
//********
lpane.add(promptLabel);
lpane.add(inputLabel);
lpane.add(promptXLabel);
lpane.add(inputXLabel);
lpane.add(promptYLabel);
lpane.add(inputYLabel);
pane1.add(lpane,BorderLayout.SOUTH);
//*********
jta=new Plot2D();
inputXmin.setText(Double.toString(jta.p1.xmin));
inputXmax.setText(Double.toString(jta.p1.xmax));
inputYmin.setText(Double.toString(jta.p1.ymin));
inputYmax.setText(Double.toString(jta.p1.ymax));
inputXntic.setText(Integer.toString(jta.p1.xntic));
inputYntic.setText(Integer.toString(jta.p1.yntic));
inputLabel.setText(jta.p1.label);
inputXLabel.setText(jta.p1.xlabel);
inputYLabel.setText(jta.p1.ylabel);
printButton=new JButton("Yazdır");
pane2.add(jta);
pane1.add(printButton,BorderLayout.SOUTH);
tabbedPane.addTab(girdi, pane3);
tabbedPane.addTab(leastcikti, pane4);
tabbedPane.addTab(PLOTPanel, pane2);
tabbedPane.addTab(KONTROLPanel, pane1);
contentPane.add( tabbedPane, BorderLayout.CENTER);
printButton.addActionListener(this);
inputXmin.addActionListener(this);
inputXmax.addActionListener(this);
inputYmin.addActionListener(this);
inputYmax.addActionListener(this);
inputLabel.addActionListener(this);
inputXLabel.addActionListener(this);
inputYLabel.addActionListener(this);
inputXntic.addActionListener(this);
inputYntic.addActionListener(this);
inputXgridon.addItemListener(this);
inputYgridon.addItemListener(this);
//en küçük kareler metodu action listener bağlantısı
inputP1_2.addActionListener(this); // nCikti
inputP1_3.addActionListener(this); // minX
inputP1_4.addActionListener(this); // maxX
inputGirdi.addActionListener(this); // Girdi dosyasi
inputCikti.addActionListener(this); // Cikti dosyasi
inputPD.addActionListener(this); // PolinomDerecesi
leastCiktisiniPlotGirdisiOlarakYazdir();
}
//===========================================================
public void itemStateChanged(ItemEvent e)
{
inputXmin.setText(Double.toString(jta.p1.xmin));
inputXmax.setText(Double.toString(jta.p1.xmax));
inputYmin.setText(Double.toString(jta.p1.ymin));
inputYmax.setText(Double.toString(jta.p1.ymax));
inputXntic.setText(Integer.toString(jta.p1.xntic));
inputYntic.setText(Integer.toString(jta.p1.yntic));
Object source=e.getItemSelectable();
if(source==inputXgridon)
{
if (e.getStateChange() == ItemEvent.DESELECTED)
{
jta.p1.xgridon=0;
}
else
{
jta.p1.xgridon=1;
}
}
else if(source==inputYgridon)
{
if (e.getStateChange() == ItemEvent.DESELECTED)
{
jta.p1.ygridon=0;
}
else
{
jta.p1.ygridon=1;
}
}
inputLabel.setText(jta.p1.label);
inputXLabel.setText(jta.p1.xlabel);
inputYLabel.setText(jta.p1.ylabel);
jta.yenidenciz();
}
//===========================================================
public void leastCiktisiniPlotGirdisiOlarakYazdir()
{
try
{
BufferedWriter fplot=new BufferedWriter(new FileWriter(plotDosyasi));
fplot.println(inputLabel.getText());
fplot.println(inputXLabel.getText());
fplot.println(inputYLabel.getText());
fplot.println(inputGirdi.getText()+" 20 0 0 0 ");
fplot.println(inputCikti.getText()+" 3 0 0 255 ");
fplot.close();
}
catch(IOException e1)
{
System.err.println("girdi cikti hatasi : Plot.txt");
}
catch(AccessControlException ace)
{
System.err.println("least ciktisini plot girdisi olarak yazdir access kontrol hatası : "+ plotDosyasi);
}
}
public void sayfalariyenile()
{
//en küçük kareler metoduna göre değişik
// sayfalardaki de§erleri yeniler
girdiDosyasi=new File(inputGirdi.getText());
Integer d1=new Integer(inputPD.getText());
PolinomDerecesi=d1.intValue();
ciktiDosyasi=new File(inputCikti.getText());
Integer d2=new Integer(inputP1_2.getText());
nCikti=d2.intValue();
Double d3=new Double(inputP1_3.getText());
minX=d3.doubleValue();
Double d4=new Double(inputP1_4.getText());
maxX=d4.doubleValue();
try{
oku();
}
catch(IOException e1)
{
System.err.println("least girdi dosyası hatası : "+ girdiDosyasi);
}
catch(AccessControlException ace)
{
System.err.println("least girdi access kontrol hatası : "+ girdiDosyasi);
}
ee=new double[PolinomDerecesi+1];
ee=Numeric.poly_least_square(xi,yi,PolinomDerecesi);
outputGirdi.setText(yazxiyi());
outputCikti.setText(yazCikti());
outputFormul.setText(yaze());
jta.yenidanPlotDatasiOku();
leastCiktisiniPlotGirdisiOlarakYazdir();
}
//===========================================================
public void actionPerformed(ActionEvent e)
{
// en kucuk kareler metodu girdileri verildiginde bu yeni action
// eventi oluşturur
if( e.getSource()==inputGirdi ||
e.getSource()==inputPD ||
e.getSource()==inputCikti ||
e.getSource()==inputP1_2 ||
e.getSource()==inputP1_3 ||
e.getSource()==inputP1_4 )
{
sayfalariyenile();
}
//grafiği yazıcıya gönderir
//********************
if (e.getSource()==printButton )
{
jta.yazdir();
}
//grafik min max değerlerini yeniden okur
Double valXmin=new Double(inputXmin.getText());
jta.p1.xmin=valXmin.doubleValue();
Double valXmax=new Double(inputXmax.getText());
jta.p1.xmax=valXmax.doubleValue();
Double valYmin=new Double(inputYmin.getText());
jta.p1.ymin=valYmin.doubleValue();
Double valYmax=new Double(inputYmax.getText());
//*****
//grafik tic sayısını yeniden okur
Integer valXntic=new Integer(inputXntic.getText());
jta.p1.xntic=valXntic.intValue();
Integer valYntic=new Integer(inputYntic.getText());
jta.p1.yntic=valYntic.intValue();
//*****
//grafik başlıklarını yeniden okur
jta.p1.ymax=valYmax.doubleValue();
jta.p1.label=inputLabel.getText();
jta.p1.xlabel=inputXLabel.getText();
jta.p1.ylabel=inputYLabel.getText();
//Plot'ı (grafik) yeniden çizer
jta.yenidenciz();
}
//===========================================================
//girdi dosyasından girdi değerlerini okur.
public void oku() throws IOException
{
//dosyadan bilgileri okur xi,yi boyutlu matrisine y�kler.
try{
fin=new BufferedReader(new FileReader(girdiDosyasi));
}
catch(FileNotFoundException e2)
{
System.err.println(e2.toString());
}
nGirdi=0;
Vector v1,v2;
v1=new Vector(1);
v2=new Vector(1);
String X1,X2;
boolean EOF=false;
while(!EOF)
{
try{
X1=Text.readString(fin);
X2=Text.readString(fin);
v1.addElement(X1);
v2.addElement(X2);
nGirdi++;
}catch(EOFException e)
{
fin.close();
EOF=true;
}
}
xi=new double[nGirdi];
yi=new double[nGirdi];
Enumeration enum1=v1.elements();
Enumeration enum2=v2.elements();
StringBuffer b=new StringBuffer();
//for(int i=0;i<n;i++)
int i=0;
while(enum1.hasMoreElements())
{
X1=(String)enum1.nextElement();
Double d1=new Double(X1);
X2=(String)enum2.nextElement();
Double d2=new Double(X2);
xi[i]=d1.doubleValue();
yi[i]=d2.doubleValue();
i++;
}
} //oku metodu sonu
//===========================================================
//girdi alanına girdi de§erlerini yazar
public String yazxiyi()
{
String b=new String();
for(int i=0;i<xi.length;i++)
{
b=b+" "+xi[i]+"\t "+yi[i]+"\n";
}
return b;
}
//========================================================
// Polinom katsayılarını yazar
public String yaze()
{
String b="Polinom katsayilari : ";
b=b+"f(x) = a[0]+a[1]*x+a[2]*x^2+....a[n]*x^n \n";
for(int i=0;i<=PolinomDerecesi;i++)
{
b=b+"a["+i+"] = "+ee[i]+"\n";
}
return b;
}
//=========================================================
// cikti dosyasina kayit yapar ve cikti alanina sonuclar� yazar
public String yazCikti()
{
String b=new String();
try{
ObjectOutputStream fout=new ObjectOutputStream(new FileOutputStream(ciktiDosyasi));
b="cikti dosyasi : \n";
double x,y;
double nn=((maxX-minX)/(nCikti-1));
for(int i=0;i<nCikti;i++)
{
x=minX+nn*i;
y=Numeric.f_least_square(ee,x);
b=b+" "+x+"\t "+y+"\n";
}
fout.writeObject(b);
try{
fout.close();
}
catch(IOException io)
{
System.exit(1);
}
}
catch(IOException e1)
{
System.err.println("input output error");
}
catch(AccessControlException ece)
{
System.err.println("yaz Cikti dosyasi"+ciktiDosyasi+ " excess kontrol");
}
return b;
}
//===========================================================
public static void main( String argv[] ) throws IOException {
leastSW frame = new leastSW( "En küçük kareler metodu eğri uydurma" );
frame.addWindowListener(new BasicWindowMonitor());
frame.setSize( 600, 600 );
frame.setVisible(true);
}
}
12001.JPG
12002.JPG
12003.JPG
12004.JPG
Şekil 12.1-4 en küçük kareler metoduyla polinaom eğri uyduran leastSW.java sınıfı Jframe çıktılarının çeşitli sayfalarının görünümü.
Örnek çoktıdaki girdi dosyaları :
başlık
x ekseni
y ekseni
2
in.txt 20 0 0 0
out.txt 3 0 0 255
in.txt
1.0 0.9999999999999432
2.0 3.9999999999997726
3.0 8.999999999999488
4.0 15.99999999999909
5.0 24.99999999999858
6.0 35.999999999997954
7.0 48.999999999997215
8.0 63.99999999999636
9.0 80.9999999999954
10.0 99.99999999999432
out.txt
cikti dosyasi :
1.0 0.9999999999999432
2.0 3.9999999999997726
3.0 8.999999999999488
4.0 15.99999999999909
5.0 24.99999999999858
6.0 35.999999999997954
7.0 48.999999999997215
8.0 63.99999999999636
9.0 80.9999999999954
10.0 99.99999999999432
Diferansiyel
Denklemler :
Problem
12.11 diferansiyel denklem testi, RK4 metodu
dif1.java
//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;
import Numeric;
import complex;
import Matrix;
class fm1 extends f_xi
{
//multivariable function
double func(double x[],int x_ref)
{
//x[0] is x, x[1] is y
//this is a representation of equation : dy/dx=y
//solution of this equation is e(x)
if(x_ref==0) return x[1];
else return 0;
}
}
class dif1
{
public static void main(String args[]) throws IOException
{
fm1 b3=new fm1();
double x[];
x=new double[1];
x[0]=1;
//RK4 differential equation to solve differential equation dy/dx=y
//with limits 0 to 1 exact solution is e=2.7182818
//boundary condition x=0 y=1 is given
//final value x=1 to be calculated by the method
System.out.println(Matrix.toStringT(Numeric.RK4(b3,0.0,1.0,1.0,50)));
}
}
Çözüm seti :
1.0
1.02020134
1.0408107741377959
1.0618365464618167
1.0832870675613175
1.1051709179307265
1.1274968514019572
1.1502737986460576
1.1735108707455981
1.197217362839226
1.2214027578398445
1.2460767302279048
1.271249149921327
1.2969300862235986
1.3231298118516308
1.3498588070449815
1.3771277637580916
1.4049475899372086
1.4333294138837107
1.4622845887055762
1.4918246968587778
1.5219615547804188
1.5527072176154666
1.5840739840389706
1.6160744011756962
1.648721269619143
1.682027648551951
1.7160068609697496
1.7506724990105325
1.7860384293916938
1.8221187989569014
1.8589280403350215
1.896480877713363
1.9347923327275491
1.9738777304703714
2.0137527056220317
2.0544332087042223
2.095935512460547
2.1382762183658373
2.1814722632669596
2.2255409261577848
2.270499835091013
2.3163669742296302
2.363160691040814
2.4108997036351645
2.4596031082541976
2.5092903869090972
2.5599814151737794
2.611696470135386
2.6644562385053905
2.7182818248945586
Problem 12.12 : çok
boyutlu RKF45 diferansiyel denklem çözüm metodu
örneği, dif2.java
//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;
import Numeric;
import complex;
import Matrix;
class fm1 extends f_xi
{
//multivariable function
double func(double x[],int x_ref)
{
//x[0] is x, x[1] is y
if(x_ref==0) return -2.0*x[0]-x[1];
else return 0;
}
}
class dif2
{
public static void main(String args[]) throws IOException
{
fm1 b3=new fm1();
double x[];
x=new double[1];
x[0]=1;
//RK5 differential equation to solve equation dy/dx=-2x-y
//initial value is given as x=0 y=-1 at x=0.4 y will be determined
System.out.println
(Matrix.toString(Matrix.T(Numeric.RKF45(b3,0.0,0.4,-1.0,10))));
}
}
Çözüm :
0.0 -1.0
0.0050 -0.9950377826066799
0.01 -0.9901501844020889
0.015 -0.9853368332473136
0.02 -0.9805973588593618
0.025 -0.9759313928019077
0.030000000000000002-0.9713385684760812
0.035 -0.9668185211113052
0.04 -0.9623708877561767
0.045 -0.9579953072693951
0.049999999999999996-0.9536914203107341
0.05499999999999999-0.94945886933206
0.05999999999999999-0.9452972985683943
0.06499999999999999-0.9412063540290199
0.06999999999999999-0.9371856834886335
0.075 -0.9332349364785408
0.08 -0.9293537642778962
0.085 -0.9255418199049855
0.09000000000000001-0.9217987581085532
0.09500000000000001-0.9181242353591716
0.10000000000000002-0.9145179098406548
0.10500000000000002-0.9109794414415133
0.11000000000000003-0.9075084917464531
0.11500000000000003-0.9041047240279163
0.12000000000000004-0.9007678032376644
0.12500000000000003-0.8974973959984025
0.13000000000000003-0.8942931705954467
0.13500000000000004-0.891154796968432
0.14000000000000004-0.8880819467030621
0.14500000000000005-0.8850742930228999
0.15000000000000005-0.8821315107811996
0.15500000000000005-0.8792532764527788
0.16000000000000006-0.8764392681259321
0.16500000000000006-0.8736891654943832
0.17000000000000007-0.8710026498492798
0.17500000000000007-0.8683794040712255
0.18000000000000008-0.8658191126223544
0.18500000000000008-0.8633214615384427
0.19000000000000009-0.8608861384210608
0.1950000000000001-0.8585128324297647
0.2000000000000001-0.856201234274326
0.2050000000000001-0.8539510362070005
0.2100000000000001-0.8517619320148359
0.2150000000000001-0.849633617012017
0.2200000000000001-0.8475657880322502
0.22500000000000012-0.8455581434211851
0.23000000000000012-0.843610383028874
0.23500000000000013-0.8417222082022695
0.24000000000000013-0.8398933217777589
0.24500000000000013-0.8381234280737365
0.2500000000000001-0.8364122328832122
0.2550000000000001-0.8347594434664576
0.2600000000000001-0.8331647685436885
0.2650000000000001-0.8316279182877833
0.27000000000000013-0.8301486043170393
0.27500000000000013-0.8287265396879636
0.28000000000000014-0.8273614388881004
0.28500000000000014-0.8260530178288946
0.29000000000000015-0.8248009938385901
0.29500000000000015-0.8236050856551641
0.30000000000000016-0.8224650134192967
0.30500000000000016-0.8213804986673754
0.31000000000000016-0.8203512643245341
0.31500000000000017-0.8193770346977273
0.3200000000000002-0.8184575354688391
0.3250000000000002-0.8175924936878259
0.3300000000000002-0.8167816377658933
0.3350000000000002-0.8160246974687078
0.3400000000000002-0.8153214039096416
0.3450000000000002-0.8146714895430517
0.3500000000000002-0.8140746881575915
0.3550000000000002-0.8135307348695573
0.3600000000000002-0.8130393661162659
0.3650000000000002-0.812600319649468
0.3700000000000002-0.8122133345287916
0.3750000000000002-0.8118781511152201
0.3800000000000002-0.8115945110646018
0.38500000000000023-0.8113621573211925
0.39000000000000024-0.8111808341112298
0.39500000000000024-0.8110502869365394
0.4 -0.8109702625681743
Problem
12.13 diferansiyel denklem çözüm örneği, RK4 metodu,
dif3.java
//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;
import Numeric;
import complex;
import Matrix;
class fm1
extends f_xi
{
//multivariable function
double func(double x[],int x_ref)
{
//x[0] is x, x[1] is y
if(x_ref==0) return -2.0*x[0]-x[1];
else return 0;
}
}
class dif3
{
public static void main(String args[]) throws IOException
{
fm1 b3=new
fm1();
double x[];
x=new double[1];
x[0]=1;
//RK4 differential equation to solve equation dy/dx=-2x-y
System.out.println(Matrix.toStringT(Numeric.RK4(b3,0.0,0.4,-1.0,100)));
System.in.read();
}
}
Çözüm :
-1.0
-0.9960239680320001
-0.992095744511233
-0.9882151385858677
-0.9843819601659565
-0.9805960199203918
-0.9768571292738785
-0.9731651004039156
-0.9695197462377918
-0.9659208804495921
-0.9623683174572164
-0.9588618724194105
-0.955401361232808
-0.9519866005289858
-0.9486174076715289
-0.945293600753109
-0.942014998592574
-0.9387814207320486
-0.9355926874340476
-0.9324486196785992
-0.9293490391603817
-0.9262937682858702
-0.9232826301704945
-0.9203154486358103
-0.9173920482066786
-0.9145122541084599
-0.9116758922642163
-0.9088827892919272
-0.9061327725017148
-0.903425669893081
-0.900761310152156
-0.898139522648956
-0.8955601374346547
-0.8930229852388627
-0.8905278974669202
-0.888074706197199
-0.8856632441784156
-0.8832933448269559
-0.8809648422242088
-0.8786775711139124
-0.8764313668995095
-0.8742260656415141
-0.872061504054888
-0.8699375195064296
-0.8678539500121697
-0.8658106342347813
-0.8638074114809979
-0.8618441216990421
-0.8599206054760645
-0.8580367040355941
-0.8561922592349969
-0.854387113562946
-0.8526211101369017
-0.8508940927006006
-0.8492059056215563
-0.8475563938885688
-0.8459454031092447
-0.8443727795075264
-0.8428383699212321
-0.8413420217996048
-0.8398835832008722
-0.8384629027898148
-0.8370798298353451
-0.8357342142080957
-0.8344259063780173
-0.8331547574119857
-0.8319206189714199
-0.8307233433099079
-0.8295627832708429
-0.8284387922850687
-0.8273512243685349
-0.826299934119961
-0.8252847767185101
-0.8243056079214711
-0.8233622840619519
-0.8224546620465801
-0.8215825993532133
-0.8207459540286591
-0.8199445846864039
-0.8191783505043502
-0.8184471112225643
Şu ana kadarki
örneklerimizin hepsi konsol programı olarak verildi. simdi grafik (applet)
çıktısına da bir örnek olmak üzere ninci dereceden denklemin köklerini
hesaplayan kokN.java
programını verelim.
Problem 12.14 ninici dereceden polinomun köklerini bulma, kokN.java programı
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import Matrix;
public class kokN extends
Applet implements ActionListener
{
private Label prompt1,prompt2;
private TextField input;
TextArea t;
Panel YaziPaneli;
int n;
String s;
public void init()
{
setBackground(Color.lightGray);
YaziPaneli=new Panel();
YaziPaneli.setFont(new Font("Serif",Font.BOLD,12));
YaziPaneli.setLayout( new GridLayout(3,1) );
t=new TextArea(8,47);
prompt1= new Label("a[0]+a[1]*x+...+a[n]*x^n=0");
prompt2= new Label("n inci dereceden polinomun katsayilarini
giriniz : ");
input = new TextField(30);
YaziPaneli.add(prompt1);
YaziPaneli.add(prompt2);
YaziPaneli.add(input);
add(YaziPaneli);
add(t);
input.addActionListener(this);
}
public void actionPerformed(ActionEvent
e)
{
s=input.getText();
StringTokenizer token=new StringTokenizer(s);
t.setText("");
n=token.countTokens()-1;
int m=n+1;
double a[]=new double[m];
complex z[]=new complex[n];
for(int i=0;i<n;i++)
{
z[i]=new complex();
}
int j=0;
while(token.hasMoreTokens())
{
Double ax=new Double(token.nextToken());
a[j++]=ax.doubleValue();
}
z=Matrix.poly_rootsC(a);
int i=0;
t.setText(Matrix.toStringT(z));
input.setText("");
}
}
Bu progamın sonucu :
12005.JPG
Şekil 12.5 n inci dereceden polinomun köklerini bulan kokN programı, StringTokenizer ve Matrix sınıf örneği
12.5 ALIŞTIRMALAR
1. f(x)=x*x*x-3 denkleminin köklerini newton metodunu kullanarak çözünüz (3 ün küp kökünü hesaplayınız), Bunun için bir java programı yazıp çıktıyı hesaplatınız.
- f(x)=x*x*x-3 denkleminin köklerini poly_roots metodunu kullanarak çözünüz.
3. f(x)=x*x*x-3 denkleminin köklerini bisection metodunu kullanarak çözünüz (3 ün küp kökünü hesaplayınız), bunun için bir java programı yazıp çıktıyı hesaplatınız.
4. f(x)=sin(x) denkleminin 0 ila pi arasındaki integralini hesaplayınız. Bunun için bir java programı yazıp çıktıyı hesaplatınız.
5. f(x)=sin(x) denkleminin pi noktasındaki türevini alınız. Bunun için bir java programı yazıp çıktıyı hesaplatınız.
6. f(x,y)=x*x+y*y-R*R (R=sabit) denkleminin 0 ile 1 arasındaki türevini hesaplayınız.
7. df(x)/dx=-f(x) denkleminin x=1 den 2 ye kadar çözümünü yapınız, f(0)=1, Bunun için bir java programı yazıp çıktıyı hesaplatınız.
<< index
< bölüm 11
> bölüm 13
bölüm başı