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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.differentiation;
23  
24  import java.io.Serializable;
25  
26  import org.hipparchus.analysis.UnivariateFunction;
27  import org.hipparchus.analysis.UnivariateMatrixFunction;
28  import org.hipparchus.analysis.UnivariateVectorFunction;
29  import org.hipparchus.exception.LocalizedCoreFormats;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.util.FastMath;
32  import org.hipparchus.util.MathArrays;
33  
34  /** Univariate functions differentiator using finite differences.
35   * <p>
36   * This class creates some wrapper objects around regular
37   * {@link UnivariateFunction univariate functions} (or {@link
38   * UnivariateVectorFunction univariate vector functions} or {@link
39   * UnivariateMatrixFunction univariate matrix functions}). These
40   * wrapper objects compute derivatives in addition to function
41   * values.
42   * </p>
43   * <p>
44   * The wrapper objects work by calling the underlying function on
45   * a sampling grid around the current point and performing polynomial
46   * interpolation. A finite differences scheme with n points is
47   * theoretically able to compute derivatives up to order n-1, but
48   * it is generally better to have a slight margin. The step size must
49   * also be small enough in order for the polynomial approximation to
50   * be good in the current point neighborhood, but it should not be too
51   * small because numerical instability appears quickly (there are several
52   * differences of close points). Choosing the number of points and
53   * the step size is highly problem dependent.
54   * </p>
55   * <p>
56   * As an example of good and bad settings, lets consider the quintic
57   * polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}.
58   * Since it is a polynomial, finite differences with at least 6 points
59   * should theoretically recover the exact same polynomial and hence
60   * compute accurate derivatives for any order. However, due to numerical
61   * errors, we get the following results for a 7 points finite differences
62   * for abscissae in the [-10, 10] range:
63   * <ul>
64   *   <li>step size = 0.25, second order derivative error about 9.97e-10</li>
65   *   <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
66   *   <li>step size = 1.0e-6, second order derivative error about 148</li>
67   *   <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li>
68   * </ul>
69   * <p>
70   * This example shows that the small step size is really bad, even simply
71   * for second order derivative!</p>
72   *
73   */
74  public class FiniteDifferencesDifferentiator
75      implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
76                 UnivariateMatrixFunctionDifferentiator, Serializable {
77  
78      /** Serializable UID. */
79      private static final long serialVersionUID = 20120917L;
80  
81      /** Number of points to use. */
82      private final int nbPoints;
83  
84      /** Step size. */
85      private final double stepSize;
86  
87      /** Half sample span. */
88      private final double halfSampleSpan;
89  
90      /** Lower bound for independent variable. */
91      private final double tMin;
92  
93      /** Upper bound for independent variable. */
94      private final double tMax;
95  
96      /**
97       * Build a differentiator with number of points and step size when independent variable is unbounded.
98       * <p>
99       * Beware that wrong settings for the finite differences differentiator
100      * can lead to highly unstable and inaccurate results, especially for
101      * high derivation orders. Using very small step sizes is often a
102      * <em>bad</em> idea.
103      * </p>
104      * @param nbPoints number of points to use
105      * @param stepSize step size (gap between each point)
106      * @exception MathIllegalArgumentException if {@code stepsize <= 0} (note that
107      * {@link MathIllegalArgumentException} extends {@link MathIllegalArgumentException})
108      * @exception MathIllegalArgumentException {@code nbPoint <= 1}
109      */
110     public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
111         throws MathIllegalArgumentException {
112         this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
113     }
114 
115     /**
116      * Build a differentiator with number of points and step size when independent variable is bounded.
117      * <p>
118      * When the independent variable is bounded (tLower &lt; t &lt; tUpper), the sampling
119      * points used for differentiation will be adapted to ensure the constraint holds
120      * even near the boundaries. This means the sample will not be centered anymore in
121      * these cases. At an extreme case, computing derivatives exactly at the lower bound
122      * will lead the sample to be entirely on the right side of the derivation point.
123      * </p>
124      * <p>
125      * Note that the boundaries are considered to be excluded for function evaluation.
126      * </p>
127      * <p>
128      * Beware that wrong settings for the finite differences differentiator
129      * can lead to highly unstable and inaccurate results, especially for
130      * high derivation orders. Using very small step sizes is often a
131      * <em>bad</em> idea.
132      * </p>
133      * @param nbPoints number of points to use
134      * @param stepSize step size (gap between each point)
135      * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
136      * if there are no lower bounds)
137      * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
138      * if there are no upper bounds)
139      * @exception MathIllegalArgumentException if {@code stepsize <= 0} (note that
140      * {@link MathIllegalArgumentException} extends {@link MathIllegalArgumentException})
141      * @exception MathIllegalArgumentException {@code nbPoint <= 1}
142      * @exception MathIllegalArgumentException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
143      */
144     public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
145                                            final double tLower, final double tUpper)
146             throws MathIllegalArgumentException {
147 
148         if (nbPoints <= 1) {
149             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
150                                                    stepSize, 1);
151         }
152         this.nbPoints = nbPoints;
153 
154         if (stepSize <= 0) {
155             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
156                                                    stepSize, 0);
157         }
158         this.stepSize = stepSize;
159 
160         halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
161         if (2 * halfSampleSpan >= tUpper - tLower) {
162             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
163                                                    2 * halfSampleSpan, tUpper - tLower);
164         }
165         final double safety = FastMath.ulp(halfSampleSpan);
166         this.tMin = tLower + halfSampleSpan + safety;
167         this.tMax = tUpper - halfSampleSpan - safety;
168 
169     }
170 
171     /**
172      * Get the number of points to use.
173      * @return number of points to use
174      */
175     public int getNbPoints() {
176         return nbPoints;
177     }
178 
179     /**
180      * Get the step size.
181      * @return step size
182      */
183     public double getStepSize() {
184         return stepSize;
185     }
186 
187     /**
188      * Evaluate derivatives from a sample.
189      * <p>
190      * Evaluation is done using divided differences.
191      * </p>
192      * @param t evaluation abscissa value and derivatives
193      * @param t0 first sample point abscissa
194      * @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)}
195      * @param <T> the type of the field elements
196      * @return value and derivatives at {@code t}
197      * @exception MathIllegalArgumentException if the requested derivation order
198      * is larger or equal to the number of points
199      */
200     private <T extends Derivative<T>> T evaluate(final T t, final double t0, final double[] y)
201         throws MathIllegalArgumentException {
202 
203         // create divided differences diagonal arrays
204         final double[] top    = new double[nbPoints];
205         final double[] bottom = new double[nbPoints];
206 
207         for (int i = 0; i < nbPoints; ++i) {
208 
209             // update the bottom diagonal of the divided differences array
210             bottom[i] = y[i];
211             for (int j = 1; j <= i; ++j) {
212                 bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize);
213             }
214 
215             // update the top diagonal of the divided differences array
216             top[i] = bottom[0];
217 
218         }
219 
220         // evaluate interpolation polynomial (represented by top diagonal) at t
221         T interpolation = t.getField().getZero();
222         T monomial      = null;
223         for (int i = 0; i < nbPoints; ++i) {
224             if (i == 0) {
225                 // start with monomial(t) = 1
226                 monomial = t.getField().getOne();
227             } else {
228                 // monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1))
229                 final T deltaX = t.subtract(t0 + (i - 1) * stepSize);
230                 monomial = monomial.multiply(deltaX);
231             }
232             interpolation = interpolation.add(monomial.multiply(top[i]));
233         }
234 
235         return interpolation;
236 
237     }
238 
239     /** {@inheritDoc}
240      * <p>The returned object cannot compute derivatives to arbitrary orders. The
241      * value function will throw a {@link MathIllegalArgumentException} if the requested
242      * derivation order is larger or equal to the number of points.
243      * </p>
244      */
245     @Override
246     public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) {
247         return new UnivariateDifferentiableFunction() {
248 
249             /** {@inheritDoc} */
250             @Override
251             public double value(final double x) throws MathIllegalArgumentException {
252                 return function.value(x);
253             }
254 
255             /** {@inheritDoc} */
256             @Override
257             public <T extends Derivative<T>> T value(T t)
258                 throws MathIllegalArgumentException {
259 
260                 // check we can achieve the requested derivation order with the sample
261                 if (t.getOrder() >= nbPoints) {
262                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
263                                                            t.getOrder(), nbPoints);
264                 }
265 
266                 // compute sample position, trying to be centered if possible
267                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
268 
269                 // compute sample points
270                 final double[] y = new double[nbPoints];
271                 for (int i = 0; i < nbPoints; ++i) {
272                     y[i] = function.value(t0 + i * stepSize);
273                 }
274 
275                 // evaluate derivatives
276                 return evaluate(t, t0, y);
277 
278             }
279 
280         };
281     }
282 
283     /** {@inheritDoc}
284      * <p>The returned object cannot compute derivatives to arbitrary orders. The
285      * value function will throw a {@link MathIllegalArgumentException} if the requested
286      * derivation order is larger or equal to the number of points.
287      * </p>
288      */
289     @Override
290     public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) {
291         return new UnivariateDifferentiableVectorFunction() {
292 
293             /** {@inheritDoc} */
294             @Override
295             public double[]value(final double x) throws MathIllegalArgumentException {
296                 return function.value(x);
297             }
298 
299             /** {@inheritDoc} */
300             @Override
301             public <T extends Derivative<T>> T[] value(T t)
302                 throws MathIllegalArgumentException {
303 
304                 // check we can achieve the requested derivation order with the sample
305                 if (t.getOrder() >= nbPoints) {
306                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
307                                                            t.getOrder(), nbPoints);
308                 }
309 
310                 // compute sample position, trying to be centered if possible
311                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
312 
313                 // compute sample points
314                 double[][] y = null;
315                 for (int i = 0; i < nbPoints; ++i) {
316                     final double[] v = function.value(t0 + i * stepSize);
317                     if (i == 0) {
318                         y = new double[v.length][nbPoints];
319                     }
320                     for (int j = 0; j < v.length; ++j) {
321                         y[j][i] = v[j];
322                     }
323                 }
324 
325                 // evaluate derivatives
326                 final T[] value = MathArrays.buildArray(t.getField(), y.length);
327                 for (int j = 0; j < value.length; ++j) {
328                     value[j] = evaluate(t, t0, y[j]);
329                 }
330 
331                 return value;
332 
333             }
334 
335         };
336     }
337 
338     /** {@inheritDoc}
339      * <p>The returned object cannot compute derivatives to arbitrary orders. The
340      * value function will throw a {@link MathIllegalArgumentException} if the requested
341      * derivation order is larger or equal to the number of points.
342      * </p>
343      */
344     @Override
345     public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) {
346         return new UnivariateDifferentiableMatrixFunction() {
347 
348             /** {@inheritDoc} */
349             @Override
350             public double[][]  value(final double x) throws MathIllegalArgumentException {
351                 return function.value(x);
352             }
353 
354             /** {@inheritDoc} */
355             @Override
356             public <T extends Derivative<T>> T[][] value(T t)
357                 throws MathIllegalArgumentException {
358 
359                 // check we can achieve the requested derivation order with the sample
360                 if (t.getOrder() >= nbPoints) {
361                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
362                                                            t.getOrder(), nbPoints);
363                 }
364 
365                 // compute sample position, trying to be centered if possible
366                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
367 
368                 // compute sample points
369                 double[][][] y = null;
370                 for (int i = 0; i < nbPoints; ++i) {
371                     final double[][] v = function.value(t0 + i * stepSize);
372                     if (i == 0) {
373                         y = new double[v.length][v[0].length][nbPoints];
374                     }
375                     for (int j = 0; j < v.length; ++j) {
376                         for (int k = 0; k < v[j].length; ++k) {
377                             y[j][k][i] = v[j][k];
378                         }
379                     }
380                 }
381 
382                 // evaluate derivatives
383                 final T[][] value = MathArrays.buildArray(t.getField(), y.length, y[0].length);
384                 for (int j = 0; j < value.length; ++j) {
385                     for (int k = 0; k < y[j].length; ++k) {
386                         value[j][k] = evaluate(t, t0, y[j][k]);
387                     }
388                 }
389 
390                 return value;
391 
392             }
393 
394         };
395     }
396 
397 }