Polynomial Fittings
Now that we have introduced concepts in Linear Algebra, it is time we apply
it to a category of numerical analysis known as Polynomial Fittings.
In the purpose of determining the behavior of a function, as data pairs
[ xi, y(xi) ] is severalfold. If we want to
approximate other values of the function at values of x not tabulated
(interpolation or extrapolation) and to estimate the integral of
f(x) and its derivative, Polynomial Fittings provides the capability
of smoothing the data. Approximating unknown values of a function is straightforward.
We will try to find a polynomial that fits a selected set of points
[xi, y(xi) ] and assume that the polynomial and
the function behave nearly the same over the interval in question. Values of the
polynomial then should be reasonable estimates of the unknown function. There are problems
with interpolating polynomials when data are not smooth, meaning there are local
irregularities. In such cases polynomial of higher order would be required to follow
the irregularities, but such polynomials, while fitting to the irregularity, deviate
widely at other regions where the function is smooth. An interpolating polynomial,
although passing through the points used in its construction, does not in general
give exactly correct values when used for interpolation, and this leads to error of
approximation.
Java Codes
Polyfit Class
This class has one constructor that instantiates a Polyfit object. The constructor
takes two input argument arrays with an integer to specify the polynomial order to
fit data into. All instance variables are private and all instance methods are
accessor. The code is descriptive as much as possible for the non-expert reader.
package jamlab;
import Jama.*;
import Jama.util.*;
public class Polyfit
{
//instance variable for upper triangular matrix
private Matrix R;
//array for polynomial coefficients
private double[] polynomialCoefficients ;
//degree of freedom
private int degreeOfFreedom;
//norm – largest singular value.
private double norm;
//y-intercept : f(x)-> x=0, f(x=0)
private double yIntercept ;
//store polynomial coefficients in matrix
private Matrix polyCoeffMatrix ;
//Contructor with arrays of data pairs , _x & _y with approximating order
public Polyfit(double[] _x , double[] _y , int order) throws Exception{
int xLength = _x.length,
yLength = _y.length;
double tempValue = 0.0;
//store array _y as a matrix
Matrix y2D = new Matrix(JElmat.convertTo2D(_y));
//transpose matrix y2D
Matrix yTranspose = y2D.transpose();
//throw Exception if data pairs are not the same size.
if(xLength != yLength){
throw new Exception(” Polyfit :- The lengths of the 2-input array parameters must be equal.”);
}
//throw Exception if only one data pair.
if(xLength < 2){
throw new Exception(" Polyfit :- There must be at least 2 data points for polynomial fitting.");
}
//throw Exception if polynomial order is negative.
if(order < 0){
throw new Exception(" Polyfit :- The polynomial fitting order cannot be less than zero.");
}
//throw Exception if polynomial order exceeds or equal the number of data-pairs.
if(order >= xLength){
throw new Exception(” Polyfit :- The polynomial order = “+order+” , must be less than the number of data points = “+xLength);
}
Matrix tempMatrix = null;
try{
//Create Vandermonde matrix
tempMatrix = JElmat.vander(_x,order+1);
}
catch(Exception ex){
ex.printStackTrace();
}
//Factorise matrix tempMatrix by QR decomposition
QRDecomposition qr = new QRDecomposition(tempMatrix);
//Initialize Orthonormal basis matrix Q from object qr
Matrix Q = qr.getQ();
//Create Upper-triangular invertible matrix
R = qr.getR();
//result = R-1 * ( QT * _yT )
Matrix result = R.inverse().times(Q.transpose().times(yTranspose));
//return a reference to the internal array transpose of the result matrix
double[][] arrayResult = result.transpose().getArray();
//polynomial coefficients is the first row of the internl array.
polynomialCoefficients = arrayResult[0];
//degree of freedom
degreeOfFreedom = yLength – (order+1);
Matrix r = yTranspose.minus(tempMatrix.times(result));
//largest singular value of matrix r
norm = r.norm2();
//Create matrix of the coefficients.
polyCoeffMatrix = new Matrix(JElmat.convertTo2D(polynomialCoefficients));
//Polynomial y-axis intercept
yIntercept = polynomialCoefficients[polynomialCoefficients.length-1];
}//end constructor
//return polynomial coefficients array
public double[] getPolynomialCoefficients(){
return polynomialCoefficients ;
}
//return upper-triangular invertible matrix
public Matrix getR(){
return R;
}
//return the degree of freedom
public int getDegreeOfFreedom(){
return degreeOfFreedom;
}
//return norm
public double getNorm(){
return norm;
}
//return y-intercept
public double getYIntercept(){
return yIntercept;
}
//return polynomial coefficients matrix
public Matrix getPolyCoeffMatrix(){
return polyCoeffMatrix ;
}
}//—– End Class Definition ——-
Polyval Class
package jamlab;
import Jama.*;
import Jama.util.*;
public class Polyval
{
//Interpolated and extrapolated values corresponding to domain _x .
private double[] Yout;
//Interpolated and extrapolated values in a matrix object .
private Matrix youtMatrix ;
//Error estimates in an array.
private double[] ErrorBounds;
//Error estimates in a matrix
private Matrix ErrorBoundsMatrix ;
//Constructor with an array of values for interpolation and extrapolation based on coefficients form polyfit object
public Polyval( double[] _x , Polyfit polyfit){
//Retrieve polynomial coefficients from polyfit object
double[] polyCoeff = polyfit.getPolynomialCoefficients();
//Length of the coefficients array
int coeffLen = polyCoeff.length;
//polyfit object degree of freedom
int degree = polyfit.getDegreeOfFreedom();
//polyfit object norm
double normr = polyfit.getNorm();
//Retrieve polyfit upper triangular matrix instance
Matrix R = polyfit.getR();
//Convert the domain _x into a matrix
Matrix _xMatrix = new Matrix(JElmat.convertTo2D(_x));
//Initialize interpolated and extrapolated output to a matrix of zeros of size 1 x length(_x)
youtMatrix = JElmat.zeros(1,_xMatrix.getColumnDimension());
Matrix pC;
//Loop to calculate estimate output extrapolated or interpolated values
for(int i=0 ; i
//youtMatrix = (youtMatrix .* _xMatrix) + pC
youtMatrix.arrayTimesEquals(_xMatrix).plusEquals(pC);
}
//Order of polynomial
int order = coeffLen;
//Declare Vandermonde matrix;
Matrix V = null;
//Declare error-bound matrix
Matrix errorDelta = null;
//Vandermonde matrix with number of columns = order
try{ V = JElmat.vander(_x,order); }
catch(Exception ex){}
//Rinv = R-1
Matrix Rinv = R.inverse();
//E = V * R-1
Matrix E = V.times(Rinv);
//Declare a matrix of ones.
Matrix matOnes = null;
//Declare a matrix for evaluation of error.
Matrix e = null;
//Horizontal line (Order = 0)
if(order == 0){
//matrix of ones with size = [E.getRowDimension(),E.getColumnDimension()]
matOnes = JElmat.ones(E.getRowDimension(),E.getColumnDimension());
//e = sqrt(matOnes + E.*E)
e = JElfun.sqrt( matOnes.plus( E.arrayTimes(E)));
}
else{ //Higher order, 1,2,3,…
//SumT = (sum((E.*E)T))T
Matrix SumT = JDatafun.sum(E.arrayTimes(E).transpose()).transpose();
//row size
int rowDim = SumT.getRowDimension();
//column size
int colDim = SumT.getColumnDimension();
//matrix of ones , same size as SumT
matOnes = JElmat.ones(rowDim,colDim);
//e = sqrt(matOnes + SumT)
e = JElfun.sqrt(matOnes.plus(SumT));
}
if(degree==1){//When degree of freedom = 1, error-bounds would be infinite.
double[][] inf = new double[1][1];
inf[0][0] = Double.POSITIVE_INFINITY ;
try {
//error matrix of size =[e.getRowDimension(),e.getColumnDimension()] with all entries is an infinity.
errorDelta = new Matrix( JElmat.repmat(inf,e.getRowDimension(),e.getColumnDimension()) );
}
catch(Exception ex) {
ex.printStackTrace();
}
}
else{ // Higher degree of freedom
errorDelta = e.times(normr/Math.sqrt(degree));
}
double[][] tempDelta = {{0.},{0.}} ; //initialize
try{
//reshape matrix errorDelta into tempDelta array with size = [1 , _x.length]
tempDelta = JElmat.reshape(errorDelta,1,_x.length).getArrayCopy() ;
}
catch(Exception ex){ ex.printStackTrace(); }
//Assign ErrorBounds first row of tempDelta
ErrorBounds = tempDelta[0] ;
//ErrorBounds in a matrix
ErrorBoundsMatrix = new Matrix(JElmat.convertTo2D(ErrorBounds));
//get copy of the output array from output matrix
double[][] yM = youtMatrix.getArrayCopy();
//Assign output to first row of the array
Yout = yM[0];
}//End constructor
//return output matrix
public Matrix getYoutMatrix(){
return youtMatrix ;
}
//return output array
public double[] getYout(){
return Yout;
}
//return output error-bounds array
public double[] getErrorBounds(){
return ErrorBounds;
}
//return output error-bounds matrix
public Matrix getErrorBoundsMatrix(){
return ErrorBoundsMatrix;
}
}// — End Class Definition —
We will conclude by running the complete program in Part 4.
Download Jama and
Jamlab and
Jamlab documentation.
References
- Java for Engineers and Scientists by Stephen J. Chapman, Prentice Hall, 1999.
- Introductory Java for Scientists and Engineers by Richard J. Davies, Addison-Wesley Pub. Co., 1999.
- Applied Numerical Analysis (Sixth Edition) by Curtis F. Gerald and Patrick O. Wheatly, Addison-Wesley Pub. Co., 1999.
- Linear Algebra and Its Applications (Second Edition), by David C. Lay, Addison-Wesley Pub. Co.
- Numerical Recipes in Fortran 77, the Art of Scientific Computing (Volume 1) by William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery, Cambridge University Press, 1997.
- Mastering MATLAB 5: A Comprehensive Tutorial and Reference by Duane Hanselman and Bruce Littlefield, Prentice Hall, 1997.
- Advanced Mathematics and Mechanics Applications using MATLAB by Louis H. Turcotte and Howard B. Wilson, CRC Press, 1998.
Related Articles
- Java as a Scientific Programming Language (Part 1): More Issues for Scientific Programming in Java
- Scientific Computing in Java (Part 2): Writing Scientific Programs in Java
- Using Java in Scientific Research: An Introduction
About the Author
Sione Palu is a Java developer at Datacom in Auckland, New Zealand, currently involved in a Web application development project. Palu graduated from the University of Auckland, New Zealand, double majoring in mathematics and computer science. He has a personal interest in applying Java and mathematics in the fields of mathematical modeling and simulations, expert systems, neural and soft computation, wavelets, digital signal processing, and control systems.