GuidesJava in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial...

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

Developer.com content and product recommendations are editorially independent. We may make money when you click on links to our partners. Learn More.

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.


<br /> package jamlab;</p> <p>import Jama.*;<br /> import Jama.util.*;</p> <p>public class Polyfit<br /> {<br /> //instance variable for upper triangular matrix<br /> private Matrix R; </p> <p> //array for polynomial coefficients<br /> private double[] polynomialCoefficients ;</p> <p> //degree of freedom<br /> private int degreeOfFreedom; </p> <p> //norm &#8211; largest singular value.<br /> private double norm; </p> <p> //y-intercept : f(x)-> x=0, f(x=0)<br /> private double yIntercept ; </p> <p> //store polynomial coefficients in matrix<br /> private Matrix polyCoeffMatrix ; </p> <p> //Contructor with arrays of data pairs , _x &#038; _y with approximating order<br /> public Polyfit(double[] _x , double[] _y , int order) throws Exception{<br /> int xLength = _x.length,<br /> yLength = _y.length;<br /> double tempValue = 0.0;</p> <p> //store array _y as a matrix<br /> Matrix y2D = new Matrix(JElmat.convertTo2D(_y)); </p> <p> //transpose matrix y2D<br /> Matrix yTranspose = y2D.transpose(); </p> <p> //throw Exception if data pairs are not the same size.<br /> if(xLength != yLength){<br /> throw new Exception(&#8221; Polyfit :- The lengths of the 2-input array parameters must be equal.&#8221;);<br /> }</p> <p> //throw Exception if only one data pair.<br /> 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){<br /> throw new Exception(&#8221; Polyfit :- The polynomial order = &#8220;+order+&#8221; , must be less than the number of data points = &#8220;+xLength);<br /> }</p> <p> Matrix tempMatrix = null;<br /> try{<br /> //Create Vandermonde matrix<br /> tempMatrix = JElmat.vander(_x,order+1);<br /> }<br /> catch(Exception ex){<br /> ex.printStackTrace();<br /> }</p> <p> //Factorise matrix tempMatrix by QR decomposition<br /> QRDecomposition qr = new QRDecomposition(tempMatrix);</p> <p> //Initialize Orthonormal basis matrix Q from object qr<br /> Matrix Q = qr.getQ();</p> <p> //Create Upper-triangular invertible matrix<br /> R = qr.getR();<br />

        //result = R-1 * ( QT * _yT )

<br /> Matrix result = R.inverse().times(Q.transpose().times(yTranspose));</p> <p> //return a reference to the internal array transpose of the result matrix<br /> double[][] arrayResult = result.transpose().getArray();</p> <p> //polynomial coefficients is the first row of the internl array.<br /> polynomialCoefficients = arrayResult[0];</p> <p> //degree of freedom<br /> degreeOfFreedom = yLength &#8211; (order+1);</p> <p> Matrix r = yTranspose.minus(tempMatrix.times(result));</p> <p> //largest singular value of matrix r<br /> norm = r.norm2();</p> <p> //Create matrix of the coefficients.<br /> polyCoeffMatrix = new Matrix(JElmat.convertTo2D(polynomialCoefficients));</p> <p> //Polynomial y-axis intercept<br /> yIntercept = polynomialCoefficients[polynomialCoefficients.length-1];</p> <p> }//end constructor</p> <p> //return polynomial coefficients array<br /> public double[] getPolynomialCoefficients(){<br /> return polynomialCoefficients ;<br /> }</p> <p> //return upper-triangular invertible matrix<br /> public Matrix getR(){<br /> return R;<br /> }</p> <p> //return the degree of freedom<br /> public int getDegreeOfFreedom(){<br /> return degreeOfFreedom;<br /> }</p> <p> //return norm<br /> public double getNorm(){<br /> return norm;<br /> }</p> <p> //return y-intercept<br /> public double getYIntercept(){<br /> return yIntercept;<br /> }</p> <p> //return polynomial coefficients matrix<br /> public Matrix getPolyCoeffMatrix(){<br /> return polyCoeffMatrix ;<br /> }</p> <p>}//&#8212;&#8211; End Class Definition &#8212;&#8212;-<br />

Polyval Class

<br /> package jamlab;</p> <p>import Jama.*;<br /> import Jama.util.*;</p> <p>public class Polyval<br /> {<br /> //Interpolated and extrapolated values corresponding to domain _x .<br /> private double[] Yout;</p> <p> //Interpolated and extrapolated values in a matrix object .<br /> private Matrix youtMatrix ;</p> <p> //Error estimates in an array.<br /> private double[] ErrorBounds;</p> <p> //Error estimates in a matrix<br /> private Matrix ErrorBoundsMatrix ;</p> <p> //Constructor with an array of values for interpolation and extrapolation based on coefficients form polyfit object<br /> public Polyval( double[] _x , Polyfit polyfit){</p> <p> //Retrieve polynomial coefficients from polyfit object<br /> double[] polyCoeff = polyfit.getPolynomialCoefficients();</p> <p> //Length of the coefficients array<br /> int coeffLen = polyCoeff.length;</p> <p> //polyfit object degree of freedom<br /> int degree = polyfit.getDegreeOfFreedom();</p> <p> //polyfit object norm<br /> double normr = polyfit.getNorm();</p> <p> //Retrieve polyfit upper triangular matrix instance<br /> Matrix R = polyfit.getR();</p> <p> //Convert the domain _x into a matrix<br /> Matrix _xMatrix = new Matrix(JElmat.convertTo2D(_x));</p> <p> //Initialize interpolated and extrapolated output to a matrix of zeros of size 1 x length(_x)<br /> youtMatrix = JElmat.zeros(1,_xMatrix.getColumnDimension());</p> <p> Matrix pC;</p> <p> //Loop to calculate estimate output extrapolated or interpolated values<br /> 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[i]);

        //youtMatrix = (youtMatrix .* _xMatrix) + pC


<br /> youtMatrix.arrayTimesEquals(_xMatrix).plusEquals(pC);<br /> }</p> <p> //Order of polynomial<br /> int order = coeffLen;</p> <p> //Declare Vandermonde matrix;<br /> Matrix V = null;</p> <p> //Declare error-bound matrix<br /> Matrix errorDelta = null;</p> <p> //Vandermonde matrix with number of columns = order<br /> try{ V = JElmat.vander(_x,order); }<br /> catch(Exception ex){}<br />

        //Rinv = R-1


<br /> Matrix Rinv = R.inverse();<br />

        //E = V * R-1


<br /> Matrix E = V.times(Rinv);</p> <p> //Declare a matrix of ones.<br /> Matrix matOnes = null;</p> <p> //Declare a matrix for evaluation of error.<br /> Matrix e = null;</p> <p> //Horizontal line (Order = 0)<br /> if(order == 0){</p> <p> //matrix of ones with size = [E.getRowDimension(),E.getColumnDimension()]<br /> matOnes = JElmat.ones(E.getRowDimension(),E.getColumnDimension());<br />

        //e = sqrt(matOnes + E.*E)


<br /> e = JElfun.sqrt( matOnes.plus( E.arrayTimes(E)));<br /> }<br /> else{ //Higher order, 1,2,3,&#8230;<br />

        //SumT = (sum((E.*E)T))T


<br /> Matrix SumT = JDatafun.sum(E.arrayTimes(E).transpose()).transpose();</p> <p> //row size<br /> int rowDim = SumT.getRowDimension();</p> <p> //column size<br /> int colDim = SumT.getColumnDimension();</p> <p> //matrix of ones , same size as SumT<br /> matOnes = JElmat.ones(rowDim,colDim);<br />

        //e = sqrt(matOnes + SumT)


<br /> e = JElfun.sqrt(matOnes.plus(SumT));<br /> }</p> <p> if(degree==1){//When degree of freedom = 1, error-bounds would be infinite.<br /> double[][] inf = new double[1][1];<br /> inf[0][0] = Double.POSITIVE_INFINITY ;<br /> try {</p> <p> //error matrix of size =[e.getRowDimension(),e.getColumnDimension()] with all entries is an infinity.<br /> errorDelta = new Matrix( JElmat.repmat(inf,e.getRowDimension(),e.getColumnDimension()) );<br /> }<br /> catch(Exception ex) {<br /> ex.printStackTrace();<br /> }<br /> }<br /> else{ // Higher degree of freedom<br /> errorDelta = e.times(normr/Math.sqrt(degree));<br /> }</p> <p> double[][] tempDelta = {{0.},{0.}} ; //initialize</p> <p> try{</p> <p> //reshape matrix errorDelta into tempDelta array with size = [1 , _x.length]<br /> tempDelta = JElmat.reshape(errorDelta,1,_x.length).getArrayCopy() ;<br /> }<br /> catch(Exception ex){ ex.printStackTrace(); }</p> <p> //Assign ErrorBounds first row of tempDelta<br /> ErrorBounds = tempDelta[0] ;</p> <p> //ErrorBounds in a matrix<br /> ErrorBoundsMatrix = new Matrix(JElmat.convertTo2D(ErrorBounds));</p> <p> //get copy of the output array from output matrix<br /> double[][] yM = youtMatrix.getArrayCopy();</p> <p> //Assign output to first row of the array<br /> Yout = yM[0];<br /> }//End constructor</p> <p> //return output matrix<br /> public Matrix getYoutMatrix(){<br /> return youtMatrix ;<br /> }</p> <p> //return output array<br /> public double[] getYout(){<br /> return Yout;<br /> }</p> <p> //return output error-bounds array<br /> public double[] getErrorBounds(){<br /> return ErrorBounds;<br /> }</p> <p> //return output error-bounds matrix<br /> public Matrix getErrorBoundsMatrix(){<br /> return ErrorBoundsMatrix;<br /> }</p> <p>}// &#8212; End Class Definition &#8212;<br />

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

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.

Get the Free Newsletter!

Subscribe to Developer Insider for top news, trends & analysis

Latest Posts

Related Stories