package jsci.maths;
import jsci.maths.matrices.*;
import jsci.maths.polynomials.RealPolynomial;
import jsci.maths.vectors.*;
/**
* The linear math library.
* This class cannot be subclassed or instantiated because all methods are static.
* @version 2.2
* @author Mark Hale
*/
public final class LinearMath extends AbstractMath {
private LinearMath() {}
// LINEAR SYSTEM
/**
* Solves the linear system Mx=v.
* @param M a double square matrix.
* @param v a double vector.
* @return the double vector x.
*/
public static AbstractDoubleVector solve(final AbstractDoubleSquareMatrix M,final AbstractDoubleVector v) {
final int n=v.dimension();
final double array[]=new double[n];
final int pivot[]=new int[n+1];
final AbstractDoubleSquareMatrix lu[]=M.luDecompose(pivot);
int i,j;
double sum;
// forward substitution
for(i=0;i=0;i--) {
sum=array[i];
for(j=i+1;j=0;i--) {
sum=array[i];
for(j=i+1;j=0;i--) {
sum=array[i];
for(j=i+1;j Math.abs(dx)) {
temp = dx / dy;
cs[1] = 1.0 / Math.sqrt( 1.0 + temp*temp );
cs[0] = temp * cs[1];
} else {
temp = dy / dx;
cs[0] = 1.0 / Math.sqrt( 1.0 + temp*temp );
cs[1] = temp * cs[0];
}
return (cs);
}
private static double[] applyPlaneRotation(double dx, double dy, double cs, double sn) {
double [] dxy = new double[2];
dxy[0] = cs * dx + sn * dy;
dxy[1] = -sn * dx + cs * dy;
return (dxy);
}
private static AbstractDoubleVector update(AbstractDoubleVector x, int k, double[][] H, double [] s, AbstractDoubleVector [] v) {
/************************************************
* Backsolve of a triangular + sub diagonal matrix:
************************************************/
for (int i = k; i >= 0; i--) {
s[i]=s[i]/H[i][i];
for (int j = i - 1; j >= 0; j--)
s[j]=s[j]-H[j][i]*s[i];
}
for (int j = 0; j <= k; j++){
x = x.add(v[j].scalarMultiply(s[j]));
}
return(x);
}
// LEAST SQUARES
/**
* Fits an nth degree polynomial to data using the method of least squares.
* @param n the degree of the polynomial (>= 0).
* @param data
* [0][] contains the x-series,
* [1][] contains the y-series.
*/
public static RealPolynomial leastSquaresFit(int n, final double data[][]) {
int i,j,nm1=n++;
double xsum,xysum,tmp;
final double mArray[][]=new double[n][n];
final double vArray[]=new double[n];
for(i=0;i=l;i--) {
f=s*offdiag[i];
b=c*offdiag[i];
if(Math.abs(f)>=Math.abs(g)) {
c=g/f;
r=Math.sqrt(c*c+1.0);
offdiag[i+1]=f*r;
s=1/r;
c*=s;
} else {
s=f/g;
r=Math.sqrt(s*s+1.0);
offdiag[i+1]=g*r;
c=1/r;
s*=c;
}
g=diag[i+1]-p;
r=(diag[i]-g)*s+2.0*c*b;
p=s*r;
diag[i+1]=g+p;
g=c*r-b;
}
diag[l]=diag[l]-p;
offdiag[l]=g;
offdiag[m]=0.0;
}
} while(m!=l);
}
}
/**
* Internal NETLIB tred1 routine.
* @author Richard Cannings
*/
private static void reduceSymmetric1_SquareToTridiagonal(final double matrix[][],final double diag[],final double offdiag[]) {
final int n=diag.length;
int i,j,k,l;
double f,g,h,hh,scale;
for(i=n-1;i>0;i--) {
l=i-1;
h=scale=0.0;
if(l>0) {
for(k=0;k<=l;k++)
scale+=Math.abs(matrix[i][k]);
if(scale==0.0)
offdiag[i]=matrix[i][l];
else {
for(k=0;k<=l;k++) {
matrix[i][k]/=scale;
h+=matrix[i][k]*matrix[i][k];
}
f=matrix[i][l];
g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h));
offdiag[i]=scale*g;
h-=f*g;
matrix[i][l]=f-g;
f=0.0;
for(j=0;j<=l;j++) {
g=0.0;
for(k=0;k<=j;k++)
g+=matrix[j][k]*matrix[i][k];
for(k=j+1;k<=l;k++)
g+=matrix[k][j]*matrix[i][k];
offdiag[j]=g/h;
f+=offdiag[j]*matrix[i][j];
}
hh=f/(h+h);
for(j=0;j<=l;j++) {
f=matrix[i][j];
offdiag[j]=g=offdiag[j]-hh*f;
for(k=0;k<=j;k++)
matrix[j][k]-=f*offdiag[k]+g*matrix[i][k];
}
}
} else
offdiag[i]=matrix[i][l];
diag[i]=h;
}
offdiag[0]=0.0;
for(i=0;i=l;i--) {
f=s*offdiag[i];
b=c*offdiag[i];
if(Math.abs(f)>=Math.abs(g)) {
c=g/f;
r=Math.sqrt(c*c+1.0);
offdiag[i+1]=f*r;
s=1/r;
c*=s;
} else {
s=f/g;
r=Math.sqrt(s*s+1.0);
offdiag[i+1]=g*r;
c=1/r;
s*=c;
}
g=diag[i+1]-p;
r=(diag[i]-g)*s+2.0*c*b;
p=s*r;
diag[i+1]=g+p;
g=c*r-b;
for(k=0;k0;i--) {
l=i-1;
h=scale=0.0;
if(l>0) {
for(k=0;k<=l;k++)
scale+=Math.abs(matrix[i][k]);
if(scale==0.0)
offdiag[i]=matrix[i][l];
else {
for(k=0;k<=l;k++) {
matrix[i][k]/=scale;
h+=matrix[i][k]*matrix[i][k];
}
f=matrix[i][l];
g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h));
offdiag[i]=scale*g;
h-=f*g;
matrix[i][l]=f-g;
f=0.0;
for(j=0;j<=l;j++) {
matrix[j][i]=matrix[i][j]/h;
g=0.0;
for(k=0;k<=j;k++)
g+=matrix[j][k]*matrix[i][k];
for(k=j+1;k<=l;k++)
g+=matrix[k][j]*matrix[i][k];
offdiag[j]=g/h;
f+=offdiag[j]*matrix[i][j];
}
hh=f/(h+h);
for(j=0;j<=l;j++) {
f=matrix[i][j];
offdiag[j]=g=offdiag[j]-hh*f;
for(k=0;k<=j;k++)
matrix[j][k]-=f*offdiag[k]+g*matrix[i][k];
}
}
} else
offdiag[i]=matrix[i][l];
diag[i]=h;
}
diag[0]=offdiag[0]=0.0;
for(i=0;i
Ads help maintain this website.