View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.analysis.differentiation;
18  
19  import java.io.Serializable;
20  
21  import org.apache.commons.math3.analysis.UnivariateFunction;
22  import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
23  import org.apache.commons.math3.analysis.UnivariateVectorFunction;
24  import org.apache.commons.math3.exception.MathIllegalArgumentException;
25  import org.apache.commons.math3.exception.NotPositiveException;
26  import org.apache.commons.math3.exception.NumberIsTooLargeException;
27  import org.apache.commons.math3.exception.NumberIsTooSmallException;
28  import org.apache.commons.math3.util.FastMath;
29  
30  /** Univariate functions differentiator using finite differences.
31   * <p>
32   * This class creates some wrapper objects around regular
33   * {@link UnivariateFunction univariate functions} (or {@link
34   * UnivariateVectorFunction univariate vector functions} or {@link
35   * UnivariateMatrixFunction univariate matrix functions}). These
36   * wrapper objects compute derivatives in addition to function
37   * value.
38   * </p>
39   * <p>
40   * The wrapper objects work by calling the underlying function on
41   * a sampling grid around the current point and performing polynomial
42   * interpolation. A finite differences scheme with n points is
43   * theoretically able to compute derivatives up to order n-1, but
44   * it is generally better to have a slight margin. The step size must
45   * also be small enough in order for the polynomial approximation to
46   * be good in the current point neighborhood, but it should not be too
47   * small because numerical instability appears quickly (there are several
48   * differences of close points). Choosing the number of points and
49   * the step size is highly problem dependent.
50   * </p>
51   * <p>
52   * As an example of good and bad settings, lets consider the quintic
53   * polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}.
54   * Since it is a polynomial, finite differences with at least 6 points
55   * should theoretically recover the exact same polynomial and hence
56   * compute accurate derivatives for any order. However, due to numerical
57   * errors, we get the following results for a 7 points finite differences
58   * for abscissae in the [-10, 10] range:
59   * <ul>
60   *   <li>step size = 0.25, second order derivative error about 9.97e-10</li>
61   *   <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
62   *   <li>step size = 1.0e-6, second order derivative error about 148</li>
63   *   <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li>
64   * </ul>
65   * This example shows that the small step size is really bad, even simply
66   * for second order derivative!
67   * </p>
68   * @since 3.1
69   */
70  public class FiniteDifferencesDifferentiator
71      implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
72                 UnivariateMatrixFunctionDifferentiator, Serializable {
73  
74      /** Serializable UID. */
75      private static final long serialVersionUID = 20120917L;
76  
77      /** Number of points to use. */
78      private final int nbPoints;
79  
80      /** Step size. */
81      private final double stepSize;
82  
83      /** Half sample span. */
84      private final double halfSampleSpan;
85  
86      /** Lower bound for independent variable. */
87      private final double tMin;
88  
89      /** Upper bound for independent variable. */
90      private final double tMax;
91  
92      /**
93       * Build a differentiator with number of points and step size when independent variable is unbounded.
94       * <p>
95       * Beware that wrong settings for the finite differences differentiator
96       * can lead to highly unstable and inaccurate results, especially for
97       * high derivation orders. Using very small step sizes is often a
98       * <em>bad</em> idea.
99       * </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 }