Michael Thomas Flanagan's Java Scientific Library

TriCubicInterpolation Class:     Tricubic Interpolation

     

Last update: 15 January 2011                                                                                                                              Main Page of Michael Thomas Flanagan's Java Scientific Library

This class contains the constructor and methods for performing an interpolation within a two dimensional array of data points, y = f(x1,x2,x3), using the following tricubic interpolation:

The coefficients cijk are calculated for each nearest neighbour eight point cube in the entered grid of data points and z1, z2 and z3 are the normalized values of x1, x2 and x3 for a unit cube.

This method requires the derivatives ∂y/∂x1, ∂y/∂x2, ∂y/∂x3, 2y/∂x1∂x2, 2y/∂x1∂x3, 2y/∂x2∂x3 and 3y/∂x1∂x2∂x3 for each entered grid point. If the values of these are known they may be entered with the y, x1 and x2 data. If the values of these derivatives are not known they are calculated by numerical differencing. Two numerical differencing methods are offered. The differencing may be performed using either the entered data points or using interpolated values close to the entered points after a tricubic spline interpolation on the entered data points.

See TriCubicSpline for a tricubic spline interpolation.

import directive: import flanagan.interpolation.TriCubicInterpolation;
Constructor public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, double[][][] dydx1, double[][][] dydx2, double[][][] dydx3, double[][][] d2ydx1dx2, double[][] d2ydx1dx3, double[][] d2ydx2dx3, double[][][] d3ydx1dx2dx3)
public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, int numerDiffOption)
Interpolate public double interpolate(double xx1, double xx2, double xx3)
Return the interpolated gradients public double[] getInterpolatedValues()
Grid point gradients public double[][][] getGridDydx1()
public double[][][] getGridDydx2()
public double[][][] getGridDydx3()
public double[][][] getGridD2ydx1dx2()
public double[][][] getGridD2ydx1dx3()
public double[][][] getGridD2ydx2dx3()
public double[][][] getGridD3ydx1dx2dx3()
Reset numerical differencing increment public static void resetDelta(double delta)
Rounding Error Options public static void noRoundingErrorCheck()
public static void potentialRoundingError(double potentialRoundingError)



CONSTRUCTOR

CONSTRUCTOR
public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, double[][][] dydx1, double[][][] dydx2, double[][][] dydx3, double[][][] d2ydx1dx2, double[][] d2ydx1dx3, double[][] d2ydx2dx3, double[][][] d3ydx1dx2dx3)
public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, int numerDiffOption)
Usage:                      TriCubicSInterpolation aa = new TriCubicInterpolation(x1, x2, x3, y, dydx1, dydx2, dydx3, d2ydx1dx2, d2ydx1dx3, d2ydx2dx3, d3ydx1dx2dx3);
Creates an instance of TriCubicInterpolation with its internal data arrays initialised to copies of the values in the x1, x2, y , dydx1, dydx2 and d2ydx1dx2 arrays where y is the tabulated function y = f(x1,x2), dydx1 is ∂y/∂x1, dydx2 is ∂y/∂x2, dydx3 is ∂y/∂x3, d2ydx1dx2 is 2y/∂x1∂x2, d2ydx1dx3 is 2y/∂x1∂x3, d2ydx2dx3 is 2y/∂x2∂x3 and d3ydx1dx2dx is 3y/∂x1∂x2∂x3. The y, dydx1, dydx2 and d2ydx1dx2 elements are ordered as {{{f(x1[0],x2[0],x3[0] . . . f(x1[0],x2[0],x3[r-1]} {{f(x1[0],x2[1],x3[0] . . . f(x1[0],x2[1],x3[r-1]} ......{f(x1[0],x2[q-1],x3[0] . . . f(x1[0],x2[q-1],x3[r-1]} } ... etc. where m is the number of x1 values and n is the number of x2 values, i.e. as in the order of the following nested for loops if you were reading in or calculating the tabulated data:
for(int i=0; i<x1.length; i++){
       for(int j=0; j<x2.length; j++){
              for(int k=0; k<x3.length; k++){
                     yValues[i][j][k]= 'read statement' or 'calculation of f(x1, x2, x3)';
              }
       }
}
and similarly for dydx1, dydx2 and d2ydx1dx2.

Usage:                      TriCubicInterpolation aa = new TriCubicInterpolation(x1, x2, x3, y, numerDiffOption);
Creates an instance of TriCubicInterpolation with its internal data arrays initialised to copies of the values in the x1, x2 and y arrays where y is the tabulated function y = f(x1,x2). The variables x1 and x2 are are one dimensional arrays of type double and y is a two dimensional array of type double. The y elements are ordered as {{{f(x1[0],x2[0],x3[0] . . . f(x1[0],x2[0],x3[r-1]} {{f(x1[0],x2[1],x3[0] . . . f(x1[0],x2[1],x3[r-1]} ......{f(x1[0],x2[q-1],x3[0] . . . f(x1[0],x2[q-1],x3[r-1]} } ... etc where m is the number of x1 values and n is the number of x2 values, i.e. as in the order of the following nested for loops if you were reading in or calculating the tabulated data:
for(int i=0; i<x1.length; i++){
       for(int j=0; j<x2.length; j++){
              for(int k=0; k<x3.length; k++){
                     yValues[i][j][k]= 'read statement' or 'calculation of f(x1, x2, x3)';
              }
       }
}
The derivatives ∂y/∂x1, ∂y/∂x2, ∂y/∂x3, 2y/∂x1∂x2, 2y/∂x1∂x3, 2y/∂x2∂x3 and 3y/∂x1∂x2∂x3 are calculated by numerical differencing. Two numerical differencing options are available:
  1. numerDiffOption = 0
    The gradients are calculated using only the supplied data points, i.e. as,

  2. numerDiffOption = 1
    The gradients are calculated using a tricubic spline interpolation, i.e. as,

    where h1, h2 and h3 are calculated as:
     h1 = delta times (maximum value of x1 - minimum value of x1)
     h2 = delta times (maximum value of x2 - minimum value of x2)
     h3 = delta times (maximum value of x3 - minimum value of x3)
    The default value of delta is 1 x 10-3.
    The values of the functions f(x1, x2, x3, h1, h2, h3) are calculated by interpolation of a tricubic spline of the entered data. If h1, h2 or h3 are greater than the mean separation of the relevant data points hx is replaced by the minimimum separation of the relevant data points.
    See below for details of how the value of delta may be reset.



METHODS

INTERPOLATION
public double interpolate(double xx1, double xx2, double xx3)
Usage:                      y1 = aa.interpolate(xx1, xx2, xx3);
Returns the interpolated value of y, y1, for given values of xx1, xx2 and xx3, using the y = f(x1,x2,x3) data entered via the constructor. This method may be called as often as required. The coefficients and derivatives needed for the interpolation are calculated and stored on instantiation so that they need not be recalculated on each call to this method. This method throws an IllegalArgumentException if xx1, xx2 or xx3 is outside the range of the values of x1[], x2[] or x3[] supplied to the constructor.



GET THE INTERPOLATED GRADIENTS
public double[] getInterpolatedValues()
Usage:                      values = aa.getInterpolatedValues();
Returns an array of eleven doubles which contains
for the last interpolation performed on calling the interpolate(xx1, xx2) method.



GET THE GRID POINT GRADIENTS
public double[][][] getGridDydx1()
public double[][][] getGridDydx2()
public double[][][] getGridDydx3()
public double[][][] getGridD2ydx1dx2()
public double[][][] getGridD2ydx1dx3()
public double[][][] getGridD2ydx2dx3()
public double[][][] getGridD3ydx1dx2dx3()
Usage:                      grad1 = aa.getGridDydx1();
This method returns the gradients, ∂y/∂x1, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad2 = aa.getGridDydx2();
This method returns the gradients, ∂y/∂x2, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad3 = aa.getGridDydx3();
This method returns the gradients, ∂y/∂x3, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad4 = aa.getGridD2ydx1dx2();
This method returns the gradients, 2y/∂x1∂x2, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad5 = aa.getGridD2ydx1dx3();
This method returns the gradients, 2y/∂x1∂x3, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad6 = aa.getGridD2ydx2dx3();
This method returns the gradients, 2y/∂x2∂x3, for each grid point, that were either entered via the constructor or calculated by the program (see below).

Usage:                      grad7 = aa.getGridD3ydx1dx2dx3();
This method returns the gradients, 3y/∂x1∂x2∂x3, for each grid point, that were either entered via the constructor or calculated by the program (see below).



RESET THE NUMERICAL DIFFERENCING FACTOR, delta
public static void resetDelta(double delta)
Usage:                      TricubicInterpolation.resetDelta();
If the gradients have not been entered via the constructer and the numerical differencing method using interpolation has been chosen the gradients are calculated as:

where h1, h2 and h3 are calculated as:
 h1 = delta times (maximum value of x1 - minimum value of x1)
 h2 = delta times (maximum value of x2 - minimum value of x2)
 h3 = delta times (maximum value of x3 - minimum value of x3)
The default value of delta is 1 x 10-3.
The values of the functions f(x1, x2, x3, h1, h2, h3) are calculated by interpolation of a tricubic spline of the entered data. If h1, h2 or h3 are greater than the mean separation of the relevant data points hx is replaced by the minimimum separation of the relevant data points.
This method allows the value of delta to be reset. If delta is to be reset this static method must be called before TricubicInterpolation is instantiated.



ROUNDING ERROR OPTIONS
public static void noRoundingErrorCheck()
Usage:                      TriCubicSpline.noRoundingErrorCheck();
Seveveral applications that call TriCubicSpline, e.g. plotting programs, may calculate an array of points to be fed to an instance of cubic spline that should lie between the limits initially supplied to the instance but which, due to rounding errors in the calculation of the array, may have extreme values that lie outside the limits by the an amount equal to the rounding error, e.g. an array that should lie between limits of 0 and 4 may run from 0 to 4.0000000000000003. The default option of the TriCubicSpline class is to check for such violations of the order of a rounding error and round the extreme value to the limit thus preventing an out of range exception being thrown. This method allows this default option to be ignored so that an out of range exception is thrown if any value lies outside he range of the limis no matter how small the violation is.

public static void potentialRoundingError(double potentialRoundingError)
Usage:                      TriCubicSpline.potentialRoundingError(potentialRoundingError);
Seveveral applications that call TriCubicSpline, e.g. plotting programs, may calculate an array of points to be fed to an instance of cubic spline that should lie between the limits initially supplied to the instance but which, due to rounding errors in the calculation of the array, may have extreme values that lie outside the limits by the an amount equal to the rounding error, e.g. an array that should lie between limits of 0 and 4 may run from 0 to 4.0000000000000003. The default option of the TriCubicSpline class is to check for such violations of the order of a rounding error and round the extreme value to the limit thus preventing an out of range exception being thrown. The default calculculation of the potential rounding error is the multiplication of 5x10-15 by 10 raised to the exponent of the value lying outside the limits. This method allows the value of 5x10-15 to be reset to the user's choice, potentialRoundingError.



EXAMPLE PROGAM

TriCubicExampleOne

This program uses calculated data and compares the interpolated y values obtained using This example program may be found on TriCubicExampleOne.java



OTHER CLASSES USED BY THIS CLASS

This class uses the following classes in this library:


This page was prepared by Dr Michael Thomas Flanagan