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    }