Back to article

Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 3

July 23, 2001

Review Part 2

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&#91;&#93; 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&#91;&#93; _x , double&#91;&#93; _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&#91;&#93;&#91;&#93; arrayResult = result.transpose().getArray(); //polynomial coefficients is the first row of the internl array. polynomialCoefficients = arrayResult&#91;0&#93;; //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&#91;polynomialCoefficients.length-1&#93;; }//end constructor //return polynomial coefficients array public double&#91;&#93; 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&#91;&#93; Yout; //Interpolated and extrapolated values in a matrix object . private Matrix youtMatrix ; //Error estimates in an array. private double&#91;&#93; 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&#91;&#93; _x , Polyfit polyfit){ //Retrieve polynomial coefficients from polyfit object double&#91;&#93; 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<coeffLen ; i++){ //Create a row-vector matrix of constant value, equals to current index of polyCoeff array. pC = new Matrix(1,_xMatrix.getColumnDimension(),polyCoeff&#91;i&#93;);
        //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 = &#91;E.getRowDimension(),E.getColumnDimension()&#93; matOnes = JElmat.ones(E.getRowDimension(),E.getColumnDimension());
        //e = sqrt(matOnes + E.*E)
e = JElfun.sqrt( 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(; } if(degree==1){//When degree of freedom = 1, error-bounds would be infinite. double&#91;&#93;&#91;&#93; inf = new double&#91;1&#93;&#91;1&#93;; inf&#91;0&#93;&#91;0&#93; = Double.POSITIVE_INFINITY ; try { //error matrix of size =&#91;e.getRowDimension(),e.getColumnDimension()&#93; 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&#91;&#93;&#91;&#93; tempDelta = {{0.},{0.}} ; //initialize try{ //reshape matrix errorDelta into tempDelta array with size = &#91;1 , _x.length&#93; tempDelta = JElmat.reshape(errorDelta,1,_x.length).getArrayCopy() ; } catch(Exception ex){ ex.printStackTrace(); } //Assign ErrorBounds first row of tempDelta ErrorBounds = tempDelta&#91;0&#93; ; //ErrorBounds in a matrix ErrorBoundsMatrix = new Matrix(JElmat.convertTo2D(ErrorBounds)); //get copy of the output array from output matrix double&#91;&#93;&#91;&#93; yM = youtMatrix.getArrayCopy(); //Assign output to first row of the array Yout = yM&#91;0&#93;; }//End constructor //return output matrix public Matrix getYoutMatrix(){ return youtMatrix ; } //return output array public double&#91;&#93; getYout(){ return Yout; } //return output error-bounds array public double&#91;&#93; 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.


  • 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

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.

Sitemap | Contact Us

Thanks for your registration, follow us on our social networks to keep up-to-date