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   * @version $Id: FiniteDifferencesDifferentiator.java 1420666 2012-12-12 13:33:20Z erans $
69   * @since 3.1
70   */
71  public class FiniteDifferencesDifferentiator
72      implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
73                 UnivariateMatrixFunctionDifferentiator, Serializable {
74  
75      /** Serializable UID. */
76      private static final long serialVersionUID = 20120917L;
77  
78      /** Number of points to use. */
79      private final int nbPoints;
80  
81      /** Step size. */
82      private final double stepSize;
83  
84      /** Half sample span. */
85      private final double halfSampleSpan;
86  
87      /** Lower bound for independent variable. */
88      private final double tMin;
89  
90      /** Upper bound for independent variable. */
91      private final double tMax;
92  
93      /**
94       * Build a differentiator with number of points and step size when independent variable is unbounded.
95       * <p>
96       * Beware that wrong settings for the finite differences differentiator
97       * can lead to highly unstable and inaccurate results, especially for
98       * high derivation orders. Using very small step sizes is often a
99       * <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 }