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 < t < 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 }