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