Java in Science: Data Interpolation and Extrapolation Using Numerical Methods of Polynomial Fittings

Java in Science: Data Interpolation and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 1
By Sione Palu

The value of a function g(x) at a specific set of points x₁,x₂,x₃,...,x_n is sometimes known, but there is no analytical expression for g(x) that enables a calculation of its value at an arbitrary point. An estimation of g(x) for an arbitrary x in a way that drawing a smooth curve through (x_i,g(x_i)) can be achieved by numerical analysis.

As an example of looking at a sales profit (in millions of dollars) P = [3.4, 5.1, 7.6, 9] for every two month period T = [2, 4, 6, 8] from a manufacturing company, and when one is trying to estimate P where T = 5(inside T-boundary) or T = 10(outside T-boundary) it would be inaccurate without data smoothing techniques of numerical analysis. If the desired x is bounded between [x₁,...,x_n] that is x is in between the largest and the smallest of the x_i's, this is called interpolation and if x is outside the boundaries then it is called extrapolation.

Polynomials

Interpolation and extrapolation techniques must model the function, within and beyond the boundaries, [x₁,...,x_n] by some functional form. The functional form that is going to be discussed in this tutorial is known as polynomial fittings. Polynomials of order n have the general form of:

y = a_nxⁿ + a_(n-1)x^(n-1) + ... + a₂x² + a₁x + a₀

n is an integern >= 0.a_n is a Real number.a_n is the coefficient of the nth term.Example of a 4th-order polynomial: y = -2x⁴ + 5x³ + 3x² - 3x + 7 , thus polynomial coefficients are: a₄ = -2, a₃ = 5, a₂ = 3, a₁ = -3, a₀ = 7

In this tutorial, we will define and derive mathematical operations mainly in matrix algebra for how to fit a polynomial of order n into a set of known or observed data pairs , [x_i,y_i]. A polynomial of any order is specified and then the coefficients are returned as an array of doubles (Java data-type).

Linear (straight-line) fitting is a polynomial of order 1 (also known as a first-order polynomial, and it has wide applications in modeling from simulations, graphics to financial forecast . Before we delve in to the Java coding , I want to introduce the concept of linear algebra and matrix operations so the general reader will understand the code with minimum difficulties.

Matrix Algebra

CAD (computer-aided design) and CAM (computer-aided manufacturing) software have revolutionized the automotive industry, for instance. Today computer graphics forms the heart, and linear algebra the soul of modern car design. Before a new car is built, engineers design and construct a mathematical car -- a wire-frame model that exists only in computer memory and on graphics display terminals. The car is also mathematically tested (to simulate impact crash-force, engine heat conduction, metal vibration-fatigue frequency, etc.) before even a prototype is built. Linear algebra is indeed the foundation of any graphics software because all the manipulations of screen images are accomplished via algebraic techniques.

If A is an m x n matrix -- that is, a matrix with m rows and n columns -- then the scalar entry in the ith row and the jth column of A is denoted by a_ij and is called the (i,j) -- entry of A. The class Matrix.java has an internal array storage as a private instance variable named A and it has an accessor method to read its value.

Jama and Jamlab

Jama is a basic linear algebra package for Java. The classes in this package enable the construction and manipulation of real, dense matrices. Jama is sufficient enough to provide functionality for routine problems and packaged in a way that is natural and understandable to non-experts. It is intended to serve as the standard matrix class for Java and will be proposed as such to the Java Grande Forum, which represents those from academia and corporations drafting proposals to Sun on Java APIs for scientific and engineering applications. Jama was developed by the National Institute of Standards and Technology (NIST) and The MathWorks Inc.. It is available for free download.

The Jama package has six classes:

• Matrix
• CholeskyDecomposition
• LUDecomposition
• QRDecomposition
• SingularValueDecomposition
• EigenvalueDecomposition

The Matrix class provides the fundamental operations of numerical linear algebra. Various constructors create matrices from two-dimensional arrays of double-precision floating point numbers. Various gets and sets provide access to submatrices and matrix elements. The basic arithmetic operations include matrix addition and multiplication, matrix norms and selected element-by-element array operations. In this tutorial , we discuss matrix operations using Jama.

Jamlab (Java Matrix Laboratory) is my own developed Java matrix package, which is based on the functionality of the popular scientific and engineering language called MatLab from MathWorks. Jamlab is my own attempt to write or translate functions written in MatLab to a Java version. This package contain classes that are subclassed from class Matrix of the package Jama. There are currently six classes in Jamlab:

• Jdatafun
• Jelfun
• Jelmat
• RowColumnIndex
• Polyfit
• Polyval

Jamlab is not yet complete -- a full version will be available for free download in the future. The reader is encouraged to download Jama and Jamlab prior to continuing on to the next part of the series.

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.

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

Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 2
By Sione Palu

Review Part 1

Sums and Scalar Multiples

Two matrices would be equal if they have the same size (that is, the same number of rows and the same number of columns) and also their corresponding entries are equal. Matrix A is not equal to matrix B, although they have the same matrix size = 3 x 2. The corresponding entries at index (1,1) are not the same for matrix A and B because A(1,1) = 3 and B(1,1) = 4, and so they are different. Note that the entry-index starts at zero.

double[][] arrayA = {{4.,0.},{5.,3.},{-2.,1.}};
//Instantiate matrix object A
Matrix A = new Matrix(arrayA);
double[][] arrayB = {{4.,0.},{5.,4.},{-2.,1.}};
//Instantiate matrix object B
Matrix B = new Matrix(arrayB);
double[][] arrayC = {{1.,7.},{7.,4.}};
//Instantiate matrix object C
Matrix C = Matrix(arrayC);
double[][] arrayD = {{1.,0.},{0.,1.}};
//Instantiate matrix object D
Matrix D = new Matrix(arrayD);

The additions (sums) and subtractions (difference) of matrices is done element-wise. The sum of A+B would result in a matrix the same size as A or B with all entries as the sum of a_ij and b_ij , for example, the resultant entry of (1,1) = 7, which is a₁₁+b₁₁ = 3+4 = 7. The sum of two different size (dimemsion) matrices is undefined such as A+C, because size(A) = 3 by 2 and size(C) = 2 by 2. In multiplying a matrix by a scaler results in a matrix with all the entries scaled by this scaler-value. 3*B is a multiplication of B by a scaler value of 3 in which all entries are now 3b_ij.

Matrix A_plus_B = A.plus(B) ;
//returned matrix has all its entries scaled by 3.
Matrix matrix_3_times_B = B.times(3.);
//This following line will throw an IllegalArgumentException because matrix dimensions do 
//not agree, matrix A is a 4 by 2 but C is a 2 by 2 size matrix. 
Matrix A_plus_C = A.plus(C);

Matrix Multiplication

Under matrix multiplication, such as A*C, the inner dimensions have to be the same for it to be defined, regardless of the outer dimensions. In this case size(A) = [3 x 2] and size(C) = [2 x 2]; therefore, matrix multiplication is defined because they have the same inner dimensions of two, [3 x 2][2 x 2]. The dimension of the result of multiplying two matrices would be the outer dimensions of the respective matrices. A*C would result in a matrix with size(A*C) = [3 x 2], that is [3 x 2][2 x 2]. Matrix multiplication results in all entries(i,j) equals the summation of element-wise multiplication of each row of A to each column of C. Multiplication of A *B is undefined, because the inner dimensions do not agree: A* B is [3 x 2]*[3 x 2]. Multiplication of matrices is not always commutative (order of operation), B *C is defined, but reversing the order as C*B is undefined but C* D = D*C. Matrix multiplication is directly in contrast with its elementary counterpart, 3 x 4 = 4 x 3 = 12. Matrix D is a special type called Identity, and this is the equivalent in elementary arithmetic of the number 1, which is the multiplication identity, 6 x 1 = 6. Identity is a square-matrix (number of rows = columns) with any dimension and all its elements at the diagonal, from top-left to bottom-right are ones and zeros everywhere. Any matrix that is multiplied to the Identity is always that matrix, examples are: A*D = A, B* D = B and C*D = C.

Matrix A_times_C = A.times(C) ;	
//The following line will throw an IllegalArgumentException because inner matrix 
//dimensions do not agree.
Matrix A_times_B = A.times(B);   

Transpose of a Matrix

Given an m x n matrix E, the transpose of E is the n x m matrix denoted by E^T, whose columns are formed from the corresponding rows of E. General matrix transposition operations:

1. (F + G)^T = F^T + G^T

Matrix result1 = F.plus(G).transpose() ;						
Matrix result2 = F.transpose().plus(G.transpose());
//Matrix result1 and result2 are equals (same value in all entries)

2. (F^T)^T = F

Matrix result3 = F.transpose().transpose() ;										
//Matrix result3 and F are equals (same value in all entries)

3. (F * G)^T = G^T * F^T

Matrix result4 = F.times(G).transpose() ;
Matrix result5 = G.transpose().times(F.transpose()) ;							
//Matrix result4 and result5 are equals (same value in all entries)

4. For any scalar r, (r * F)^T = r * F^T

Matrix result6 = F.times(r).transpose() ;
Matrix result7 = F.transpose().times(r) ;									
//Matrix result6 and result7 are equals (same value in all entries)

Inverse of a Matrix

If K is an [n x n] matrix, sometimes there is another [n x n] matrix L such that:

K * L = I and L * K = I
where I is an [n x n] Identity matrix

We say K is invertible and L is an inverse of K. The inverse of K is denoted by K^-1.

K * K^-1 = I, where I - is the Identity matrix
K^-1 * K = I

Matrix inverseMatrix = K.inverse() ;												
//Matrix inverseMatrix is an inverse of matrix K

QR Matrix Factorization

Factorizing a matrix follows a similar concept in elementary arithmetics. The number 30 has prime factors of = 2 x 3 x 5. If P is an [m x n] matrix with linearly independent columns, then it can be factored as:

P = Q * R:
where Q is an [m x m] whose columns form an orthonormal basis and R is an [m x n] upper triangular invertible matrix with positive entries on its diagonal.

double[][] arrayP = {{1.,0.},{1.,1.},{1.,1.},{1.,1.}};
Matrix P = new Matrix(arrayP);

//Create a QR factorisation object from matrix P.
QRDecomposition QR = new QRDecomposition(P);
Matrix Q = QR.getQ();
Matrix R = QR.getR(); 
//Now that we have two matrix factors for matrix P which:  Q*R = P

Elementary Matrix Function and Special Matrices

• repmat (Matrix mat, int row, int col) - static method from JElmat class of Jamlab package
  Repeat a matrix in a tiling manner. If matrix A is to be tiled [3 by 2], then A would be tiled in a 3-rows by 2-columns pattern.

  

  double[][] arrayA = {{7,5},{2,3}}; 
  Matrix A = new Matrix(arrayA);
  //Create matrix object B by tiling matrix A in a 3 by 2 pattern.
  Matrix B = JElmat.repmat(A,3,2); 
  

• vander (double[] a,int numCol) - static method from JElmat.
  Create a Vandermonde matrix in which the second to last column are the elements of array a[]. The matrix has numCol number of columns.

  

  //Create Vandermonde-matrix
  double[] arrayA = {5,-1,7,6};
  //vander is also overloaded with one input argument a double[] array, where columns = rows
  Matrix A = JElmat.vander(arrayA);
  double[] arrayB = {5,-1,7,6,-3,4};
  //Create a Vandermonde-matrix with only four columns.
  Matrix B = JElmat.vander(arrayB,4);
  

• sum (Matrix matrix) - static method from JDatafun class.
  sum all the elements of the column of a matrix and return a single row matrix.

  

  	
  double arrayA = {{19,18,16,18},{5,15,9,15},{12,9,12,4},{10,0,16,8}};
  Matrix A = new Matrix(arrayA);
  //Sum all the columns of matrix A.
  Matrix B = JDatafun.sum(A); 
  

• linespace (double leftBound, double rightBound, int nPoints) - static method from JElmat class
  Linear spaced matrix row vector. Generates nPoints between leftBound and rightBound.

  //row matrix A, starting at -8.	and ending at 12. with eleven points equally spaced.
  //internal array A = (-8    -6    -4    -2     0     2     4     6     8    10    12)
  Matrix A = JElmat.linspace(-8.,12.,11); 
  

• reshape (Matrix A, int M, int N) - static method from JElmat class.
  returns the M-by-N matrix whose elements are taken columnwise from A.

  

  double[][] arrayA = {{1.,2.,3.,4.},{5.,6.,7.,8.},{9.,10.,11.,12.}}; 
  Matrix A = new Matrix(arrayA); 
  //reshape matrix A from 3 x 4 dimension into 6 x 2
  Matrix B = JElmat.reshape(A,6,2);
  

• zeros (int m, int n) - static method from JElmat class.
  returns an m x n matrix of zeros.
  
  

• ones (int m, int n) - static method from JElmat class.
  returns an m x n matrix of ones.

  

  //Create a matrix of ones with 3 x 4 size.
  Matrix X = JElmat.ones(3,4);
  //Create a matrix of zeros with 2 x 2 size.
  Matrix Y = JElmat.zeros(2,2);
  

  We will continue with polynomial fittings and Java code in Part 3.

  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.

  Previous article: Java in Science: Data Interpolation and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 1
  Next article: Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 3
  

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

  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 [ x_i, y(x_i) ] 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 [x_i, y(x_i) ] 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<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

           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.

  Previous article: Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 2
  Next article: Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 4
  

  Java in Science: Data Inter- and Extrapolation Using Numerical Methods of Polynomial Fittings, Part 4
  By Sione Palu
  

  Review Part 3

  Running the Complete Program

  Here is the complete program. Please review the lessons of Parts 1 through 3 before using it.

            import Jama.*; import jamlab.*; public Class Main { public static void main(String[] args) throws Exception { //actual array data pairs (xVal, yVal) double[] xVal = {-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.,12.}; double[] yVal = {-4.,-7.,-2.,1.,2.,1.5,4.,5.,4.9,6.7,9.}; //the total number of polynomials to be fitted to the data pairs int MAXIMUM_ORDER = 8; //store the polynomial coefficients double[][] Coef = new double[MAXIMUM_ORDER][]; //upper error-bounds double[][] U = new double[MAXIMUM_ORDER][]; //lower error-bounds double[][] L = new double[MAXIMUM_ORDER][]; //evaluation of best fit polynomial double[][] yEval = new double[MAXIMUM_ORDER][]; //order of polynomial int ORDER = 0; //number of points for smoothing the interpolation and extrapolation of data pairs int number_of_points = 2000; //Array of Polyfit Polyfit[] polyf = new Polyfit[MAXIMUM_ORDER]; //array of Polyval Polyval[] polyv = new Polyval[MAXIMUM_ORDER]; //array of matrices for error estimation for each polynomial order Matrix[] error = new Matrix[MAXIMUM_ORDER]; //array of matrices for polynomial evaluation for each polynomial order Matrix[] yOutMat = new Matrix[MAXIMUM_ORDER]; //array of matrices for upper-bound error estimation for each polynomial order Matrix[] upper = new Matrix[MAXIMUM_ORDER]; //array of matrices for lower-bound error estimation for each polynomial order Matrix[] lower = new Matrix[MAXIMUM_ORDER]; try{ for(int i=0; i < MAXIMUM_ORDER ; i++){ ORDER = i ; //instantiate polyfit objects with different order, starting at order=0 polyf[i] = new Polyfit(xVal,yVal,ORDER); //array of polynomial coefficients at order = ORDER Coef[ORDER] = polyf.getPolyCoeffMatrix().getColumnPackedCopy(); } } catch(Exception ex) { throw new Exception(ex.getMessage());} Matrix xm = JElmat.linspace(-8,22,number_of_points); double[] xEval = new JElmat(xm.getArrayCopy()).matrixArrayRow(0); for(int k=0 ; k < MAXIMUM_ORDER ; k++){ polyv[k] = new Polyval(xEval,polyf[k]); //matrix estimate of errorBounds error[k] = polyv[k].getErrorBoundsMatrix(); //matrix evaluation of polynomial fitting yOutMat[k] = polyv[k].getYoutMatrix(); //matrix estimate of upper-limit errorBounds upper[k] = yOutMat[k].plus(error[k]); //matrix estimate of lower-limit errorBounds lower[k] = yOutMat.minus(error[k]); //array estimate of upper-limit errorBounds U[k] = new JElmat( upper[k].getArrayCopy()).matrixArrayRow(0); //array estimate of lower-limit errorBounds L[k] = new JElmat( lower[k].getArrayCopy()).matrixArrayRow(0); //array evaluation of polynomial fitting yEval[k] = polyv[k].getYout(); } //Readers who use Java2D graphics for plotting, can implement the following loop // to plot the output evaluation and error-bounds (upper and lower). for(int j=0 ; j < MAXIMUM_ORDER ; j++) { // (xVal,yVal) -> plot actual data pairs // (xEval,yEval[j]) -> plot evaluation data pairs of best fit polynomial with order = j // (xEval,L[j]) -> plot lower error-bounds for best fit polynomial of order = j // (xEval,U[j]) -> plot upper error-bounds for best fit polynomial of order = j } }//-- end main -- }//--- End Class Definition --- 
          

  Interpolations

  Intermediate values between data-pairs can be estimated for any polynomial order by evaluating the specific polynomial in the domain boundry with a large number of points. The data pairs used in this tutorial were :
     double[] xVal = {-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.,12.};
     double[] yVal = {-4.,-7.,-2.,1.,2.,1.5,4.,5.,4.9,6.7,9.};
  There were two thousand points between the minimum , -8. and the maximum, 12, to make the estimation as smooth as possible.
  

  The following graphs show different plottings for different order of best fit polynomials.
      blue line - actual data pairs
      red line - best fit polynomial
      black line - upper bound error estimate
      yellow line - lower bound error estimate
  

  Y = 1.918182
  
  
  Y = 0.684091*X + 0.55
  
  
  Y = 0.010111*X² + 0.724534*X + 0.913986
  
  
  Y = 0.000874*X³ - 0.015355*X² + 0.672786*X + 1.031469
  
  
  Y = 0.000754*X⁴ - 0.005157*X³ - 0.072654*X² + 0.950233*X + 1.61049
  
  
  Y = 0.000131*X⁵ + 0.002061*X⁴ + 0.006171*X³ - 0.161537*X² + 0.739704*X + 2.27972
  
  
  Y = 0.000033*X⁶ - 0.000525*X⁵ - 0.000927*X⁴ + 0.040596*X³ - 0.092734*X² + 0.172263*X + 2.004319
  
  
  Y = 0.000004*X⁷ + 0.000088*X⁶ - 0.000177*X⁵ - 0.006612*X⁴ + 0.033296*X³ + 0.049072*X² + 0.205033*X + 1.506746
  

  Extrapolations

  Extending the estimation of the best fit polynomials to extrapolate outside the data pairs domain gives wildly osscilated error-bounds. Error-bounds increase as the estimation moves away from the minimum = -8 or maximum = 12 point of the data pairs domain. Error-bounds would have reasonable values around the vicinity of the minimum and the maximum.

  The following graphs plot the extrapolation beyond the data pairs domain extending from minimum = -8 to maximum = 22:

  
  
  
  
  
  
  

  Applications

  Polynomial fittings is one of the numerical techniques that engineers and scientists use most often in estimation, simulation, and even model prediction. There are other numerical methods for estimation, such as splines for example. The choice depends on the domain of interest for the engineer. A NASA scientist might prefer the cubic spline as an estimator, while a financial engineering economist in the Tokyo stock market prefer to use a polynomial fitting of ORDER = 1, (straight line) for future outcome predictions.

  Final Word

  We have introduced the concepts of matrix algebra in this tutorial and how it is applied in numerical techniques of polynomial fittings. Linear algebra is a cornerstone of today's scientific software. There are other branches of engineering mathematics, such as discrete mathematics, differential equations, multi-variable calculus, statistics and probabilities, complex variables, and operations research, etc., which are also the foundations of scientific software. I will try to visit some of these in future tutorials.

  Book References

  • Java for Engineers and Scientists by Stephen J. Chapman, publisher : PRENTICE HALL
  • Introductory Java for Scientists and Engineers by Richard J. Davies, publisher : ADDISON-WESLEY
  • Applied Numerical Analysis (Sixth Edition), by Curtis F. Gerald and Patrick O. Wheatly, publisher : ADDISON-WESLEY
  • Linear Algebra and its Applications (Second Edition), by David C. Lay, publisher : ADDISON-WESLEY
  • 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, publisher : CAMBRIDGE UNIVERSITY PRESS
  • Mastering MATLAB 5, by Duane Hanselman and Bruce Littlefield, publisher : PRENTICE HALL
  • Advanced Mathematics and Mechanics Applications using MATLAB by Howard B. Wilson , publisher : CRC Press

  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.

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