    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.differentiation;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.analysis.UnivariateFunction;
    022    import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
    023    import org.apache.commons.math3.analysis.UnivariateVectorFunction;
    024    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    025    import org.apache.commons.math3.exception.NotPositiveException;
    026    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    027    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    028    import org.apache.commons.math3.util.FastMath;
    029
    030    /** Univariate functions differentiator using finite differences.
    031     * <p>
    032     * This class creates some wrapper objects around regular
    033     * {@link UnivariateFunction univariate functions} (or {@link
    034     * UnivariateVectorFunction univariate vector functions} or {@link
    035     * UnivariateMatrixFunction univariate matrix functions}). These
    036     * wrapper objects compute derivatives in addition to function
    037     * value.
    038     * </p>
    039     * <p>
    040     * The wrapper objects work by calling the underlying function on
    041     * a sampling grid around the current point and performing polynomial
    042     * interpolation. A finite differences scheme with n points is
    043     * theoretically able to compute derivatives up to order n-1, but
    044     * it is generally better to have a slight margin. The step size must
    045     * also be small enough in order for the polynomial approximation to
    046     * be good in the current point neighborhood, but it should not be too
    047     * small because numerical instability appears quickly (there are several
    048     * differences of close points). Choosing the number of points and
    049     * the step size is highly problem dependent.
    050     * </p>
    051     * <p>
    052     * As an example of good and bad settings, lets consider the quintic
    053     * polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}.
    054     * Since it is a polynomial, finite differences with at least 6 points
    055     * should theoretically recover the exact same polynomial and hence
    056     * compute accurate derivatives for any order. However, due to numerical
    057     * errors, we get the following results for a 7 points finite differences
    058     * for abscissae in the [-10, 10] range:
    059     * <ul>
    060     *   <li>step size = 0.25, second order derivative error about 9.97e-10</li>
    061     *   <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
    062     *   <li>step size = 1.0e-6, second order derivative error about 148</li>
    063     *   <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li>
    064     * </ul>
    065     * This example shows that the small step size is really bad, even simply
    066     * for second order derivative!
    067     * </p>
    068     * @version $Id: FiniteDifferencesDifferentiator.java 1420666 2012-12-12 13:33:20Z erans $
    069     * @since 3.1
    070     */
    071    public class FiniteDifferencesDifferentiator
    072        implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
    073                   UnivariateMatrixFunctionDifferentiator, Serializable {
    074
    075        /** Serializable UID. */
    076        private static final long serialVersionUID = 20120917L;
    077
    078        /** Number of points to use. */
    079        private final int nbPoints;
    080
    081        /** Step size. */
    082        private final double stepSize;
    083
    084        /** Half sample span. */
    085        private final double halfSampleSpan;
    086
    087        /** Lower bound for independent variable. */
    088        private final double tMin;
    089
    090        /** Upper bound for independent variable. */
    091        private final double tMax;
    092
    093        /**
    094         * Build a differentiator with number of points and step size when independent variable is unbounded.
    095         * <p>
    096         * Beware that wrong settings for the finite differences differentiator
    097         * can lead to highly unstable and inaccurate results, especially for
    098         * high derivation orders. Using very small step sizes is often a
    099         * <em>bad</em> idea.
    100         * </p>
    101         * @param nbPoints number of points to use
    102         * @param stepSize step size (gap between each point)
    103         * @exception NotPositiveException if {@code stepsize <= 0} (note that
    104         * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
    105         * @exception NumberIsTooSmallException {@code nbPoint <= 1}
    106         */
    107        public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
    108            throws NotPositiveException, NumberIsTooSmallException {
    109            this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
    110        }
    111
    112        /**
    113         * Build a differentiator with number of points and step size when independent variable is bounded.
    114         * <p>
    115         * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
    116         * points used for differentiation will be adapted to ensure the constraint holds
    117         * even near the boundaries. This means the sample will not be centered anymore in
    118         * these cases. At an extreme case, computing derivatives exactly at the lower bound
    119         * will lead the sample to be entirely on the right side of the derivation point.
    120         * </p>
    121         * <p>
    122         * Note that the boundaries are considered to be excluded for function evaluation.
    123         * </p>
    124         * <p>
    125         * Beware that wrong settings for the finite differences differentiator
    126         * can lead to highly unstable and inaccurate results, especially for
    127         * high derivation orders. Using very small step sizes is often a
    128         * <em>bad</em> idea.
    129         * </p>
    130         * @param nbPoints number of points to use
    131         * @param stepSize step size (gap between each point)
    132         * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
    133         * if there are no lower bounds)
    134         * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
    135         * if there are no upper bounds)
    136         * @exception NotPositiveException if {@code stepsize <= 0} (note that
    137         * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
    138         * @exception NumberIsTooSmallException {@code nbPoint <= 1}
    139         * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
    140         */
    141        public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
    142                                               final double tLower, final double tUpper)
    143                throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
    144
    145            if (nbPoints <= 1) {
    146                throw new NumberIsTooSmallException(stepSize, 1, false);
    147            }
    148            this.nbPoints = nbPoints;
    149
    150            if (stepSize <= 0) {
    151                throw new NotPositiveException(stepSize);
    152            }
    153            this.stepSize = stepSize;
    154
    155            halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
    156            if (2 * halfSampleSpan >= tUpper - tLower) {
    157                throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
    158            }
    159            final double safety = FastMath.ulp(halfSampleSpan);
    160            this.tMin = tLower + halfSampleSpan + safety;
    161            this.tMax = tUpper - halfSampleSpan - safety;
    162
    163        }
    164
    165        /**
    166         * Get the number of points to use.
    167         * @return number of points to use
    168         */
    169        public int getNbPoints() {
    170            return nbPoints;
    171        }
    172
    173        /**
    174         * Get the step size.
    175         * @return step size
    176         */
    177        public double getStepSize() {
    178            return stepSize;
    179        }
    180
    181        /**
    182         * Evaluate derivatives from a sample.
    183         * <p>
    184         * Evaluation is done using divided differences.
    185         * </p>
    186         * @param t evaluation abscissa value and derivatives
    187         * @param t0 first sample point abscissa
    188         * @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)}
    189         * @return value and derivatives at {@code t}
    190         * @exception NumberIsTooLargeException if the requested derivation order
    191         * is larger or equal to the number of points
    192         */
    193        private DerivativeStructure evaluate(final DerivativeStructure t, final double t0,
    194                                             final double[] y)
    195            throws NumberIsTooLargeException {
    196
    197            // create divided differences diagonal arrays
    198            final double[] top    = new double[nbPoints];
    199            final double[] bottom = new double[nbPoints];
    200
    201            for (int i = 0; i < nbPoints; ++i) {
    202
    203                // update the bottom diagonal of the divided differences array
    204                bottom[i] = y[i];
    205                for (int j = 1; j <= i; ++j) {
    206                    bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize);
    207                }
    208
    209                // update the top diagonal of the divided differences array
    210                top[i] = bottom[0];
    211
    212            }
    213
    214            // evaluate interpolation polynomial (represented by top diagonal) at t
    215            final int order            = t.getOrder();
    216            final int parameters       = t.getFreeParameters();
    217            final double[] derivatives = t.getAllDerivatives();
    218            final double dt0           = t.getValue() - t0;
    219            DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0);
    220            DerivativeStructure monomial = null;
    221            for (int i = 0; i < nbPoints; ++i) {
    222                if (i == 0) {
    223                    // start with monomial(t) = 1
    224                    monomial = new DerivativeStructure(parameters, order, 1.0);
    225                } else {
    226                    // monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1))
    227                    derivatives[0] = dt0 - (i - 1) * stepSize;
    228                    final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
    229                    monomial = monomial.multiply(deltaX);
    230                }
    231                interpolation = interpolation.add(monomial.multiply(top[i]));
    232            }
    233
    234            return interpolation;
    235
    236        }
    237
    238        /** {@inheritDoc}
    239         * <p>The returned object cannot compute derivatives to arbitrary orders. The
    240         * value function will throw a {@link NumberIsTooLargeException} if the requested
    241         * derivation order is larger or equal to the number of points.
    242         * </p>
    243         */
    244        public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) {
    245            return new UnivariateDifferentiableFunction() {
    246
    247                /** {@inheritDoc} */
    248                public double value(final double x) throws MathIllegalArgumentException {
    249                    return function.value(x);
    250                }
    251
    252                /** {@inheritDoc} */
    253                public DerivativeStructure value(final DerivativeStructure t)
    254                    throws MathIllegalArgumentException {
    255
    256                    // check we can achieve the requested derivation order with the sample
    257                    if (t.getOrder() >= nbPoints) {
    258                        throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
    259                    }
    260
    261                    // compute sample position, trying to be centered if possible
    262                    final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
    263
    264                    // compute sample points
    265                    final double[] y = new double[nbPoints];
    266                    for (int i = 0; i < nbPoints; ++i) {
    267                        y[i] = function.value(t0 + i * stepSize);
    268                    }
    269
    270                    // evaluate derivatives
    271                    return evaluate(t, t0, y);
    272
    273                }
    274
    275            };
    276        }
    277
    278        /** {@inheritDoc}
    279         * <p>The returned object cannot compute derivatives to arbitrary orders. The
    280         * value function will throw a {@link NumberIsTooLargeException} if the requested
    281         * derivation order is larger or equal to the number of points.
    282         * </p>
    283         */
    284        public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) {
    285            return new UnivariateDifferentiableVectorFunction() {
    286
    287                /** {@inheritDoc} */
    288                public double[]value(final double x) throws MathIllegalArgumentException {
    289                    return function.value(x);
    290                }
    291
    292                /** {@inheritDoc} */
    293                public DerivativeStructure[] value(final DerivativeStructure t)
    294                    throws MathIllegalArgumentException {
    295
    296                    // check we can achieve the requested derivation order with the sample
    297                    if (t.getOrder() >= nbPoints) {
    298                        throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
    299                    }
    300
    301                    // compute sample position, trying to be centered if possible
    302                    final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
    303
    304                    // compute sample points
    305                    double[][] y = null;
    306                    for (int i = 0; i < nbPoints; ++i) {
    307                        final double[] v = function.value(t0 + i * stepSize);
    308                        if (i == 0) {
    309                            y = new double[v.length][nbPoints];
    310                        }
    311                        for (int j = 0; j < v.length; ++j) {
    312                            y[j][i] = v[j];
    313                        }
    314                    }
    315
    316                    // evaluate derivatives
    317                    final DerivativeStructure[] value = new DerivativeStructure[y.length];
    318                    for (int j = 0; j < value.length; ++j) {
    319                        value[j] = evaluate(t, t0, y[j]);
    320                    }
    321
    322                    return value;
    323
    324                }
    325
    326            };
    327        }
    328
    329        /** {@inheritDoc}
    330         * <p>The returned object cannot compute derivatives to arbitrary orders. The
    331         * value function will throw a {@link NumberIsTooLargeException} if the requested
    332         * derivation order is larger or equal to the number of points.
    333         * </p>
    334         */
    335        public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) {
    336            return new UnivariateDifferentiableMatrixFunction() {
    337
    338                /** {@inheritDoc} */
    339                public double[][]  value(final double x) throws MathIllegalArgumentException {
    340                    return function.value(x);
    341                }
    342
    343                /** {@inheritDoc} */
    344                public DerivativeStructure[][]  value(final DerivativeStructure t)
    345                    throws MathIllegalArgumentException {
    346
    347                    // check we can achieve the requested derivation order with the sample
    348                    if (t.getOrder() >= nbPoints) {
    349                        throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
    350                    }
    351
    352                    // compute sample position, trying to be centered if possible
    353                    final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
    354
    355                    // compute sample points
    356                    double[][][] y = null;
    357                    for (int i = 0; i < nbPoints; ++i) {
    358                        final double[][] v = function.value(t0 + i * stepSize);
    359                        if (i == 0) {
    360                            y = new double[v.length][v[0].length][nbPoints];
    361                        }
    362                        for (int j = 0; j < v.length; ++j) {
    363                            for (int k = 0; k < v[j].length; ++k) {
    364                                y[j][k][i] = v[j][k];
    365                            }
    366                        }
    367                    }
    368
    369                    // evaluate derivatives
    370                    final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length];
    371                    for (int j = 0; j < value.length; ++j) {
    372                        for (int k = 0; k < y[j].length; ++k) {
    373                            value[j][k] = evaluate(t, t0, y[j][k]);
    374                        }
    375                    }
    376
    377                    return value;
    378
    379                }
    380
    381            };
    382        }
    383
    384    }