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 */
017package org.apache.commons.math3.analysis.differentiation;
018
019import java.io.Serializable;
020
021import org.apache.commons.math3.analysis.UnivariateFunction;
022import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
023import org.apache.commons.math3.analysis.UnivariateVectorFunction;
024import org.apache.commons.math3.exception.MathIllegalArgumentException;
025import org.apache.commons.math3.exception.NotPositiveException;
026import org.apache.commons.math3.exception.NumberIsTooLargeException;
027import org.apache.commons.math3.exception.NumberIsTooSmallException;
028import 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 * @since 3.1
069 */
070public class FiniteDifferencesDifferentiator
071    implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
072               UnivariateMatrixFunctionDifferentiator, Serializable {
073
074    /** Serializable UID. */
075    private static final long serialVersionUID = 20120917L;
076
077    /** Number of points to use. */
078    private final int nbPoints;
079
080    /** Step size. */
081    private final double stepSize;
082
083    /** Half sample span. */
084    private final double halfSampleSpan;
085
086    /** Lower bound for independent variable. */
087    private final double tMin;
088
089    /** Upper bound for independent variable. */
090    private final double tMax;
091
092    /**
093     * Build a differentiator with number of points and step size when independent variable is unbounded.
094     * <p>
095     * Beware that wrong settings for the finite differences differentiator
096     * can lead to highly unstable and inaccurate results, especially for
097     * high derivation orders. Using very small step sizes is often a
098     * <em>bad</em> idea.
099     * </p>
100     * @param nbPoints number of points to use
101     * @param stepSize step size (gap between each point)
102     * @exception NotPositiveException if {@code stepsize <= 0} (note that
103     * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
104     * @exception NumberIsTooSmallException {@code nbPoint <= 1}
105     */
106    public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
107        throws NotPositiveException, NumberIsTooSmallException {
108        this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
109    }
110
111    /**
112     * Build a differentiator with number of points and step size when independent variable is bounded.
113     * <p>
114     * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
115     * points used for differentiation will be adapted to ensure the constraint holds
116     * even near the boundaries. This means the sample will not be centered anymore in
117     * these cases. At an extreme case, computing derivatives exactly at the lower bound
118     * will lead the sample to be entirely on the right side of the derivation point.
119     * </p>
120     * <p>
121     * Note that the boundaries are considered to be excluded for function evaluation.
122     * </p>
123     * <p>
124     * Beware that wrong settings for the finite differences differentiator
125     * can lead to highly unstable and inaccurate results, especially for
126     * high derivation orders. Using very small step sizes is often a
127     * <em>bad</em> idea.
128     * </p>
129     * @param nbPoints number of points to use
130     * @param stepSize step size (gap between each point)
131     * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
132     * if there are no lower bounds)
133     * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
134     * if there are no upper bounds)
135     * @exception NotPositiveException if {@code stepsize <= 0} (note that
136     * {@link NotPositiveException} extends {@link NumberIsTooSmallException})
137     * @exception NumberIsTooSmallException {@code nbPoint <= 1}
138     * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
139     */
140    public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
141                                           final double tLower, final double tUpper)
142            throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
143
144        if (nbPoints <= 1) {
145            throw new NumberIsTooSmallException(stepSize, 1, false);
146        }
147        this.nbPoints = nbPoints;
148
149        if (stepSize <= 0) {
150            throw new NotPositiveException(stepSize);
151        }
152        this.stepSize = stepSize;
153
154        halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
155        if (2 * halfSampleSpan >= tUpper - tLower) {
156            throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
157        }
158        final double safety = FastMath.ulp(halfSampleSpan);
159        this.tMin = tLower + halfSampleSpan + safety;
160        this.tMax = tUpper - halfSampleSpan - safety;
161
162    }
163
164    /**
165     * Get the number of points to use.
166     * @return number of points to use
167     */
168    public int getNbPoints() {
169        return nbPoints;
170    }
171
172    /**
173     * Get the step size.
174     * @return step size
175     */
176    public double getStepSize() {
177        return stepSize;
178    }
179
180    /**
181     * Evaluate derivatives from a sample.
182     * <p>
183     * Evaluation is done using divided differences.
184     * </p>
185     * @param t evaluation abscissa value and derivatives
186     * @param t0 first sample point abscissa
187     * @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)}
188     * @return value and derivatives at {@code t}
189     * @exception NumberIsTooLargeException if the requested derivation order
190     * is larger or equal to the number of points
191     */
192    private DerivativeStructure evaluate(final DerivativeStructure t, final double t0,
193                                         final double[] y)
194        throws NumberIsTooLargeException {
195
196        // create divided differences diagonal arrays
197        final double[] top    = new double[nbPoints];
198        final double[] bottom = new double[nbPoints];
199
200        for (int i = 0; i < nbPoints; ++i) {
201
202            // update the bottom diagonal of the divided differences array
203            bottom[i] = y[i];
204            for (int j = 1; j <= i; ++j) {
205                bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize);
206            }
207
208            // update the top diagonal of the divided differences array
209            top[i] = bottom[0];
210
211        }
212
213        // evaluate interpolation polynomial (represented by top diagonal) at t
214        final int order            = t.getOrder();
215        final int parameters       = t.getFreeParameters();
216        final double[] derivatives = t.getAllDerivatives();
217        final double dt0           = t.getValue() - t0;
218        DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0);
219        DerivativeStructure monomial = null;
220        for (int i = 0; i < nbPoints; ++i) {
221            if (i == 0) {
222                // start with monomial(t) = 1
223                monomial = new DerivativeStructure(parameters, order, 1.0);
224            } else {
225                // monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1))
226                derivatives[0] = dt0 - (i - 1) * stepSize;
227                final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
228                monomial = monomial.multiply(deltaX);
229            }
230            interpolation = interpolation.add(monomial.multiply(top[i]));
231        }
232
233        return interpolation;
234
235    }
236
237    /** {@inheritDoc}
238     * <p>The returned object cannot compute derivatives to arbitrary orders. The
239     * value function will throw a {@link NumberIsTooLargeException} if the requested
240     * derivation order is larger or equal to the number of points.
241     * </p>
242     */
243    public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) {
244        return new UnivariateDifferentiableFunction() {
245
246            /** {@inheritDoc} */
247            public double value(final double x) throws MathIllegalArgumentException {
248                return function.value(x);
249            }
250
251            /** {@inheritDoc} */
252            public DerivativeStructure value(final DerivativeStructure t)
253                throws MathIllegalArgumentException {
254
255                // check we can achieve the requested derivation order with the sample
256                if (t.getOrder() >= nbPoints) {
257                    throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
258                }
259
260                // compute sample position, trying to be centered if possible
261                final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
262
263                // compute sample points
264                final double[] y = new double[nbPoints];
265                for (int i = 0; i < nbPoints; ++i) {
266                    y[i] = function.value(t0 + i * stepSize);
267                }
268
269                // evaluate derivatives
270                return evaluate(t, t0, y);
271
272            }
273
274        };
275    }
276
277    /** {@inheritDoc}
278     * <p>The returned object cannot compute derivatives to arbitrary orders. The
279     * value function will throw a {@link NumberIsTooLargeException} if the requested
280     * derivation order is larger or equal to the number of points.
281     * </p>
282     */
283    public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) {
284        return new UnivariateDifferentiableVectorFunction() {
285
286            /** {@inheritDoc} */
287            public double[]value(final double x) throws MathIllegalArgumentException {
288                return function.value(x);
289            }
290
291            /** {@inheritDoc} */
292            public DerivativeStructure[] value(final DerivativeStructure t)
293                throws MathIllegalArgumentException {
294
295                // check we can achieve the requested derivation order with the sample
296                if (t.getOrder() >= nbPoints) {
297                    throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
298                }
299
300                // compute sample position, trying to be centered if possible
301                final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
302
303                // compute sample points
304                double[][] y = null;
305                for (int i = 0; i < nbPoints; ++i) {
306                    final double[] v = function.value(t0 + i * stepSize);
307                    if (i == 0) {
308                        y = new double[v.length][nbPoints];
309                    }
310                    for (int j = 0; j < v.length; ++j) {
311                        y[j][i] = v[j];
312                    }
313                }
314
315                // evaluate derivatives
316                final DerivativeStructure[] value = new DerivativeStructure[y.length];
317                for (int j = 0; j < value.length; ++j) {
318                    value[j] = evaluate(t, t0, y[j]);
319                }
320
321                return value;
322
323            }
324
325        };
326    }
327
328    /** {@inheritDoc}
329     * <p>The returned object cannot compute derivatives to arbitrary orders. The
330     * value function will throw a {@link NumberIsTooLargeException} if the requested
331     * derivation order is larger or equal to the number of points.
332     * </p>
333     */
334    public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) {
335        return new UnivariateDifferentiableMatrixFunction() {
336
337            /** {@inheritDoc} */
338            public double[][]  value(final double x) throws MathIllegalArgumentException {
339                return function.value(x);
340            }
341
342            /** {@inheritDoc} */
343            public DerivativeStructure[][]  value(final DerivativeStructure t)
344                throws MathIllegalArgumentException {
345
346                // check we can achieve the requested derivation order with the sample
347                if (t.getOrder() >= nbPoints) {
348                    throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
349                }
350
351                // compute sample position, trying to be centered if possible
352                final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
353
354                // compute sample points
355                double[][][] y = null;
356                for (int i = 0; i < nbPoints; ++i) {
357                    final double[][] v = function.value(t0 + i * stepSize);
358                    if (i == 0) {
359                        y = new double[v.length][v[0].length][nbPoints];
360                    }
361                    for (int j = 0; j < v.length; ++j) {
362                        for (int k = 0; k < v[j].length; ++k) {
363                            y[j][k][i] = v[j][k];
364                        }
365                    }
366                }
367
368                // evaluate derivatives
369                final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length];
370                for (int j = 0; j < value.length; ++j) {
371                    for (int k = 0; k < y[j].length; ++k) {
372                        value[j][k] = evaluate(t, t0, y[j][k]);
373                    }
374                }
375
376                return value;
377
378            }
379
380        };
381    }
382
383}