diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java | 384 |
1 files changed, 384 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java b/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java new file mode 100644 index 0000000..c2f1002 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java @@ -0,0 +1,384 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.differentiation; + +import java.io.Serializable; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.UnivariateMatrixFunction; +import org.apache.commons.math3.analysis.UnivariateVectorFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.FastMath; + +/** Univariate functions differentiator using finite differences. + * <p> + * This class creates some wrapper objects around regular + * {@link UnivariateFunction univariate functions} (or {@link + * UnivariateVectorFunction univariate vector functions} or {@link + * UnivariateMatrixFunction univariate matrix functions}). These + * wrapper objects compute derivatives in addition to function + * values. + * </p> + * <p> + * The wrapper objects work by calling the underlying function on + * a sampling grid around the current point and performing polynomial + * interpolation. A finite differences scheme with n points is + * theoretically able to compute derivatives up to order n-1, but + * it is generally better to have a slight margin. The step size must + * also be small enough in order for the polynomial approximation to + * be good in the current point neighborhood, but it should not be too + * small because numerical instability appears quickly (there are several + * differences of close points). Choosing the number of points and + * the step size is highly problem dependent. + * </p> + * <p> + * As an example of good and bad settings, lets consider the quintic + * polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}. + * Since it is a polynomial, finite differences with at least 6 points + * should theoretically recover the exact same polynomial and hence + * compute accurate derivatives for any order. However, due to numerical + * errors, we get the following results for a 7 points finite differences + * for abscissae in the [-10, 10] range: + * <ul> + * <li>step size = 0.25, second order derivative error about 9.97e-10</li> + * <li>step size = 0.25, fourth order derivative error about 5.43e-8</li> + * <li>step size = 1.0e-6, second order derivative error about 148</li> + * <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li> + * </ul> + * <p> + * This example shows that the small step size is really bad, even simply + * for second order derivative!</p> + * + * @since 3.1 + */ +public class FiniteDifferencesDifferentiator + implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator, + UnivariateMatrixFunctionDifferentiator, Serializable { + + /** Serializable UID. */ + private static final long serialVersionUID = 20120917L; + + /** Number of points to use. */ + private final int nbPoints; + + /** Step size. */ + private final double stepSize; + + /** Half sample span. */ + private final double halfSampleSpan; + + /** Lower bound for independent variable. */ + private final double tMin; + + /** Upper bound for independent variable. */ + private final double tMax; + + /** + * Build a differentiator with number of points and step size when independent variable is unbounded. + * <p> + * Beware that wrong settings for the finite differences differentiator + * can lead to highly unstable and inaccurate results, especially for + * high derivation orders. Using very small step sizes is often a + * <em>bad</em> idea. + * </p> + * @param nbPoints number of points to use + * @param stepSize step size (gap between each point) + * @exception NotPositiveException if {@code stepsize <= 0} (note that + * {@link NotPositiveException} extends {@link NumberIsTooSmallException}) + * @exception NumberIsTooSmallException {@code nbPoint <= 1} + */ + public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize) + throws NotPositiveException, NumberIsTooSmallException { + this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY); + } + + /** + * Build a differentiator with number of points and step size when independent variable is bounded. + * <p> + * When the independent variable is bounded (tLower < t < tUpper), the sampling + * points used for differentiation will be adapted to ensure the constraint holds + * even near the boundaries. This means the sample will not be centered anymore in + * these cases. At an extreme case, computing derivatives exactly at the lower bound + * will lead the sample to be entirely on the right side of the derivation point. + * </p> + * <p> + * Note that the boundaries are considered to be excluded for function evaluation. + * </p> + * <p> + * Beware that wrong settings for the finite differences differentiator + * can lead to highly unstable and inaccurate results, especially for + * high derivation orders. Using very small step sizes is often a + * <em>bad</em> idea. + * </p> + * @param nbPoints number of points to use + * @param stepSize step size (gap between each point) + * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY} + * if there are no lower bounds) + * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY} + * if there are no upper bounds) + * @exception NotPositiveException if {@code stepsize <= 0} (note that + * {@link NotPositiveException} extends {@link NumberIsTooSmallException}) + * @exception NumberIsTooSmallException {@code nbPoint <= 1} + * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower} + */ + public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize, + final double tLower, final double tUpper) + throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + + if (nbPoints <= 1) { + throw new NumberIsTooSmallException(stepSize, 1, false); + } + this.nbPoints = nbPoints; + + if (stepSize <= 0) { + throw new NotPositiveException(stepSize); + } + this.stepSize = stepSize; + + halfSampleSpan = 0.5 * stepSize * (nbPoints - 1); + if (2 * halfSampleSpan >= tUpper - tLower) { + throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false); + } + final double safety = FastMath.ulp(halfSampleSpan); + this.tMin = tLower + halfSampleSpan + safety; + this.tMax = tUpper - halfSampleSpan - safety; + + } + + /** + * Get the number of points to use. + * @return number of points to use + */ + public int getNbPoints() { + return nbPoints; + } + + /** + * Get the step size. + * @return step size + */ + public double getStepSize() { + return stepSize; + } + + /** + * Evaluate derivatives from a sample. + * <p> + * Evaluation is done using divided differences. + * </p> + * @param t evaluation abscissa value and derivatives + * @param t0 first sample point abscissa + * @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)} + * @return value and derivatives at {@code t} + * @exception NumberIsTooLargeException if the requested derivation order + * is larger or equal to the number of points + */ + private DerivativeStructure evaluate(final DerivativeStructure t, final double t0, + final double[] y) + throws NumberIsTooLargeException { + + // create divided differences diagonal arrays + final double[] top = new double[nbPoints]; + final double[] bottom = new double[nbPoints]; + + for (int i = 0; i < nbPoints; ++i) { + + // update the bottom diagonal of the divided differences array + bottom[i] = y[i]; + for (int j = 1; j <= i; ++j) { + bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize); + } + + // update the top diagonal of the divided differences array + top[i] = bottom[0]; + + } + + // evaluate interpolation polynomial (represented by top diagonal) at t + final int order = t.getOrder(); + final int parameters = t.getFreeParameters(); + final double[] derivatives = t.getAllDerivatives(); + final double dt0 = t.getValue() - t0; + DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0); + DerivativeStructure monomial = null; + for (int i = 0; i < nbPoints; ++i) { + if (i == 0) { + // start with monomial(t) = 1 + monomial = new DerivativeStructure(parameters, order, 1.0); + } else { + // monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1)) + derivatives[0] = dt0 - (i - 1) * stepSize; + final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives); + monomial = monomial.multiply(deltaX); + } + interpolation = interpolation.add(monomial.multiply(top[i])); + } + + return interpolation; + + } + + /** {@inheritDoc} + * <p>The returned object cannot compute derivatives to arbitrary orders. The + * value function will throw a {@link NumberIsTooLargeException} if the requested + * derivation order is larger or equal to the number of points. + * </p> + */ + public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) { + return new UnivariateDifferentiableFunction() { + + /** {@inheritDoc} */ + public double value(final double x) throws MathIllegalArgumentException { + return function.value(x); + } + + /** {@inheritDoc} */ + public DerivativeStructure value(final DerivativeStructure t) + throws MathIllegalArgumentException { + + // check we can achieve the requested derivation order with the sample + if (t.getOrder() >= nbPoints) { + throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); + } + + // compute sample position, trying to be centered if possible + final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; + + // compute sample points + final double[] y = new double[nbPoints]; + for (int i = 0; i < nbPoints; ++i) { + y[i] = function.value(t0 + i * stepSize); + } + + // evaluate derivatives + return evaluate(t, t0, y); + + } + + }; + } + + /** {@inheritDoc} + * <p>The returned object cannot compute derivatives to arbitrary orders. The + * value function will throw a {@link NumberIsTooLargeException} if the requested + * derivation order is larger or equal to the number of points. + * </p> + */ + public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) { + return new UnivariateDifferentiableVectorFunction() { + + /** {@inheritDoc} */ + public double[]value(final double x) throws MathIllegalArgumentException { + return function.value(x); + } + + /** {@inheritDoc} */ + public DerivativeStructure[] value(final DerivativeStructure t) + throws MathIllegalArgumentException { + + // check we can achieve the requested derivation order with the sample + if (t.getOrder() >= nbPoints) { + throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); + } + + // compute sample position, trying to be centered if possible + final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; + + // compute sample points + double[][] y = null; + for (int i = 0; i < nbPoints; ++i) { + final double[] v = function.value(t0 + i * stepSize); + if (i == 0) { + y = new double[v.length][nbPoints]; + } + for (int j = 0; j < v.length; ++j) { + y[j][i] = v[j]; + } + } + + // evaluate derivatives + final DerivativeStructure[] value = new DerivativeStructure[y.length]; + for (int j = 0; j < value.length; ++j) { + value[j] = evaluate(t, t0, y[j]); + } + + return value; + + } + + }; + } + + /** {@inheritDoc} + * <p>The returned object cannot compute derivatives to arbitrary orders. The + * value function will throw a {@link NumberIsTooLargeException} if the requested + * derivation order is larger or equal to the number of points. + * </p> + */ + public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) { + return new UnivariateDifferentiableMatrixFunction() { + + /** {@inheritDoc} */ + public double[][] value(final double x) throws MathIllegalArgumentException { + return function.value(x); + } + + /** {@inheritDoc} */ + public DerivativeStructure[][] value(final DerivativeStructure t) + throws MathIllegalArgumentException { + + // check we can achieve the requested derivation order with the sample + if (t.getOrder() >= nbPoints) { + throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); + } + + // compute sample position, trying to be centered if possible + final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; + + // compute sample points + double[][][] y = null; + for (int i = 0; i < nbPoints; ++i) { + final double[][] v = function.value(t0 + i * stepSize); + if (i == 0) { + y = new double[v.length][v[0].length][nbPoints]; + } + for (int j = 0; j < v.length; ++j) { + for (int k = 0; k < v[j].length; ++k) { + y[j][k][i] = v[j][k]; + } + } + } + + // evaluate derivatives + final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length]; + for (int j = 0; j < value.length; ++j) { + for (int k = 0; k < y[j].length; ++k) { + value[j][k] = evaluate(t, t0, y[j][k]); + } + } + + return value; + + } + + }; + } + +} |