4.1 Overview

The analysis package is the parent package for algorithms dealing with real-valued functions of one real variable. It contains dedicated sub-packages providing numerical root-finding, integration, and interpolation. It also contains a polynomials sub-package that considers polynomials with real coefficients as differentiable real functions.

Functions interfaces are intended to be implemented by user code to represent their domain problems. The algorithms provided by the library will then operate on these function to find their roots, or integrate them, or ... Functions can be multivariate or univariate, real vectorial or matrix valued, and they can be differentiable or not.

Possible future additions may include numerical differentiation.

4.2 Error handling

For user-defined functions, when the method encounters an error during evaluation, users must use their own unchecked exceptions. The following example shows the recommended way to do that, using root solving as the example (the same construct should be used for ODE integrators or for optimizations).

private static class LocalException extends RuntimeException {

   // the x value that caused the problem
   private final double x;

   public LocalException(double x) {
     this.x = x;
   }

   public double getX() {
     return x;
   }

 }

 private static class MyFunction implements UnivariateFunction {
   public double value(double x) {
     double y = hugeFormula(x);
     if (somethingBadHappens) {
       throw new LocalException(x);
     }
     return y;
   }
 }

 public void compute() {
   try {
     solver.solve(maxEval, new MyFunction(a, b, c), min, max);
   } catch (LocalException le) {
     // retrieve the x value
   }
 }
 

As shown in this example the exception is really something local to user code and there is a guarantee Apache Commons Math will not mess with it. The user is safe.

4.3 Root-finding

UnivariateSolver, DifferentiableUnivariateSolver and PolynomialSolver provide means to find roots of univariate real-valued functions, differentiable univariate real-valued functions, and polynomial functions respectively. A root is the value where the function takes the value 0. Commons-Math includes implementations of the several root-finding algorithms:

Some algorithms require that the initial search interval brackets the root (i.e. the function values at interval end points have opposite signs). Some algorithms preserve bracketing throughout computation and allow user to specify which side of the convergence interval to select as the root. It is also possible to force a side selection after a root has been found even for algorithms that do not provide this feature by themselves. This is useful for example in sequential search, for which a new search interval is started after a root has been found in order to find the next root. In this case, user must select a side to ensure his loop is not stuck on one root and always return the same solution without making any progress.

There are numerous non-obvious traps and pitfalls in root finding. First, the usual disclaimers due to the way real world computers calculate values apply. If the computation of the function provides numerical instabilities, for example due to bit cancellation, the root finding algorithms may behave badly and fail to converge or even return bogus values. There will not necessarily be an indication that the computed root is way off the true value. Secondly, the root finding problem itself may be inherently ill-conditioned. There is a "domain of indeterminacy", the interval for which the function has near zero absolute values around the true root, which may be large. Even worse, small problems like roundoff error may cause the function value to "numerically oscillate" between negative and positive values. This may again result in roots way off the true value, without indication. There is not much a generic algorithm can do if ill-conditioned problems are met. A way around this is to transform the problem in order to get a better conditioned function. Proper selection of a root-finding algorithm and its configuration parameters requires knowledge of the analytical properties of the function under analysis and numerical analysis techniques. Users are encouraged to consult a numerical analysis text (or a numerical analyst) when selecting and configuring a solver.

In order to use the root-finding features, first a solver object must be created by calling its constructor, often providing relative and absolute accuracy. Using a solver object, roots of functions are easily found using the solve methods. These methods takes a maximum iteration count maxEval, a function f, and either two domain values, min and max, or a startValue as parameters. If the maximal number of iterations count is exceeded, non-convergence is assumed and a ConvergenceException exception is thrown. A suggested value is 100, which should be plenty, given that a bisection algorithm can't get any more accurate after 52 iterations because of the number of mantissa bits in a double precision floating point number. If a number of ill-conditioned problems is to be solved, this number can be decreased in order to avoid wasting time. Bracketed solvers also take an allowed solution enum parameter to specify which side of the final convergence interval should be selected as the root. It can be ANY_SIDE, LEFT_SIDE, RIGHT_SIDE, BELOW_SIDE or ABOVE_SIDE. Left and right are used to specify the root along the function parameter axis while below and above refer to the function value axis. The solve methods compute a value c such that:

• f(c) = 0.0 (see "function value accuracy")
• min <= c <= max (except for the secant method, which may find a solution outside the interval)

Typical usage:

UnivariateFunction function = // some user defined function object
final double relativeAccuracy = 1.0e-12;
final double absoluteAccuracy = 1.0e-8;
final int    maxOrder         = 5;
UnivariateSolver solver   = new BracketingNthOrderBrentSolver(relativeAccuracy, absoluteAccuracy, maxOrder);
double c = solver.solve(100, function, 1.0, 5.0, AllowedSolution.LEFT_SIDE);

Force bracketing, by refining a base solution found by a non-bracketing solver:

UnivariateFunction function = // some user defined function object
final double relativeAccuracy = 1.0e-12;
final double absoluteAccuracy = 1.0e-8;
UnivariateSolver nonBracketing = new BrentSolver(relativeAccuracy, absoluteAccuracy);
double baseRoot = nonBracketing.solve(100, function, 1.0, 5.0);
double c = UnivariateSolverUtils.forceSide(100, function,
                                           new PegasusSolver(relativeAccuracy, absoluteAccuracy),
                                           baseRoot, 1.0, 5.0, AllowedSolution.LEFT_SIDE);

The BrentSolver uses the Brent-Dekker algorithm which is fast and robust. If there are multiple roots in the interval, or there is a large domain of indeterminacy, the algorithm will converge to a random root in the interval without indication that there are problems. Interestingly, the examined text book implementations all disagree in details of the convergence criteria. Also each implementation had problems for one of the test cases, so the expressions had to be fudged further. Don't expect to get exactly the same root values as for other implementations of this algorithm.

The BracketingNthOrderBrentSolver uses an extension of the Brent-Dekker algorithm which uses inverse n^th order polynomial interpolation instead of inverse quadratic interpolation, and which allows selection of the side of the convergence interval for result bracketing. This is now the recommended algorithm for most users since it has the largest order, doesn't require derivatives, has guaranteed convergence and allows result bracket selection.

The SecantSolver uses a straightforward secant algorithm which does not bracket the search and therefore does not guarantee convergence. It may be faster than Brent on some well-behaved functions.

The RegulaFalsiSolver is variation of secant preserving bracketing, but then it may be slow, as one end point of the search interval will become fixed after and only the other end point will converge to the root, hence resulting in a search interval size that does not decrease to zero.

The IllinoisSolver and PegasusSolver are well-known variations of regula falsi that fix the problem of stuck end points by slightly weighting one endpoint to balance the interval at next iteration. Pegasus is often faster than Illinois. Pegasus may be the algorithm of choice for selecting a specific side of the convergence interval.

The BisectionSolver is included for completeness and for establishing a fall back in cases of emergency. The algorithm is simple, most likely bug free and guaranteed to converge even in very adverse circumstances which might cause other algorithms to malfunction. The drawback is of course that it is also guaranteed to be slow.

The UnivariateSolver interface exposes many properties to control the convergence of a solver. The accuracy properties are set at solver instance creation and cannot be changed afterwards, there are only getters to retriveve their values, no setters are available.

4.4 Interpolation

A UnivariateInterpolator is used to find a univariate real-valued function f which for a given set of ordered pairs (x_i,y_i) yields f(x_i)=y_i to the best accuracy possible. The result is provided as an object implementing the UnivariateFunction interface. It can therefore be evaluated at any point, including point not belonging to the original set. Currently, only an interpolator for generating natural cubic splines and a polynomial interpolator are available. There is no interpolator factory, mainly because the interpolation algorithm is more determined by the kind of the interpolated function rather than the set of points to interpolate. There aren't currently any accuracy controls either, as interpolation accuracy is in general determined by the algorithm.

Typical usage:

double x[] = { 0.0, 1.0, 2.0 };
double y[] = { 1.0, -1.0, 2.0);
UnivariateInterpolator interpolator = new SplineInterpolator();
UnivariateFunction function = interpolator.interpolate(x, y);
double interpolationX = 0.5;
double interpolatedY = function.evaluate(x);
System.out println("f(" + interpolationX + ") = " + interpolatedY);

A natural cubic spline is a function consisting of a polynomial of third degree for each subinterval determined by the x-coordinates of the interpolated points. A function interpolating N value pairs consists of N-1 polynomials. The function is continuous, smooth and can be differentiated twice. The second derivative is continuous but not smooth. The x values passed to the interpolator must be ordered in ascending order. It is not valid to evaluate the function for values outside the range x₀..x_N.

The polynomial function returned by the Neville's algorithm is a single polynomial guaranteed to pass exactly through the interpolation points. The degree of the polynomial is the number of points minus 1 (i.e. the interpolation polynomial for a three points set will be a quadratic polynomial). Despite the fact the interpolating polynomials is a perfect approximation of a function at interpolation points, it may be a loose approximation between the points. Due to Runge's phenomenom the error can get worse as the degree of the polynomial increases, so adding more points does not always lead to a better interpolation.

Loess (or Lowess) interpolation is a robust interpolation useful for smoothing univariate scaterplots. It has been described by William Cleveland in his 1979 seminal paper Robust Locally Weighted Regression and Smoothing Scatterplots. This kind of interpolation is computationally intensive but robust.

Microsphere interpolation is a robust multidimensional interpolation algorithm. It has been described in William Dudziak's MS thesis.

A BivariateGridInterpolator is used to find a bivariate real-valued function f which for a given set of tuples (x_i,y_j,f_ij) yields f(x_i,y_j)=f_ij to the best accuracy possible. The result is provided as an object implementing the BivariateFunction interface. It can therefore be evaluated at any point, including a point not belonging to the original set. The arrays x_i and y_j must be sorted in increasing order in order to define a two-dimensional grid.

In bicubic interpolation, the interpolation function is a 3rd-degree polynomial of two variables. The coefficients are computed from the function values sampled on a grid, as well as the values of the partial derivatives of the function at those grid points. From two-dimensional data sampled on a grid, the BicubicSplineInterpolator computes a bicubic interpolating function. Prior to computing an interpolating function, the SmoothingPolynomialBicubicSplineInterpolator class performs smoothing of the data by computing the polynomial that best fits each of the one-dimensional curves along each of the coordinate axes.

A TrivariateGridInterpolator is used to find a trivariate real-valued function f which for a given set of tuples (x_i,y_j,z_k, f_ijk) yields f(x_i,y_j,z_k)=f_ijk to the best accuracy possible. The result is provided as an object implementing the TrivariateFunction interface. It can therefore be evaluated at any point, including a point not belonging to the original set. The arrays x_i, y_j and z_k must be sorted in increasing order in order to define a three-dimensional grid.

In tricubic interpolation, the interpolation function is a 3rd-degree polynomial of three variables. The coefficients are computed from the function values sampled on a grid, as well as the values of the partial derivatives of the function at those grid points. From three-dimensional data sampled on a grid, the TricubicSplineInterpolator computes a tricubic interpolating function.

4.5 Integration

A UnivariateIntegrator provides the means to numerically integrate univariate real-valued functions. Commons-Math includes implementations of the following integration algorithms:

• Romberg's method
• Simpson's method
• trapezoid method
• Legendre-Gauss method

4.6 Polynomials

The org.apache.commons.math3.analysis.polynomials package provides real coefficients polynomials.

The PolynomialFunction class is the most general one, using traditional coefficients arrays. The PolynomialsUtils utility class provides static factory methods to build Chebyshev, Hermite, Jacobi, Laguerre and Legendre polynomials. Coefficients are computed using exact fractions so these factory methods can build polynomials up to any degree.