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 */ 017package org.apache.commons.math4.legacy.analysis.differentiation; 018 019import org.apache.commons.math4.legacy.analysis.UnivariateFunction; 020import org.apache.commons.math4.legacy.analysis.UnivariateMatrixFunction; 021import org.apache.commons.math4.legacy.analysis.UnivariateVectorFunction; 022import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException; 023import org.apache.commons.math4.legacy.exception.NotPositiveException; 024import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException; 025import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 026import org.apache.commons.math4.core.jdkmath.JdkMath; 027 028/** Univariate functions differentiator using finite differences. 029 * <p> 030 * This class creates some wrapper objects around regular 031 * {@link UnivariateFunction univariate functions} (or {@link 032 * UnivariateVectorFunction univariate vector functions} or {@link 033 * UnivariateMatrixFunction univariate matrix functions}). These 034 * wrapper objects compute derivatives in addition to function 035 * values. 036 * </p> 037 * <p> 038 * The wrapper objects work by calling the underlying function on 039 * a sampling grid around the current point and performing polynomial 040 * interpolation. A finite differences scheme with n points is 041 * theoretically able to compute derivatives up to order n-1, but 042 * it is generally better to have a slight margin. The step size must 043 * also be small enough in order for the polynomial approximation to 044 * be good in the current point neighborhood, but it should not be too 045 * small because numerical instability appears quickly (there are several 046 * differences of close points). Choosing the number of points and 047 * the step size is highly problem dependent. 048 * </p> 049 * <p> 050 * As an example of good and bad settings, lets consider the quintic 051 * polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}. 052 * Since it is a polynomial, finite differences with at least 6 points 053 * should theoretically recover the exact same polynomial and hence 054 * compute accurate derivatives for any order. However, due to numerical 055 * errors, we get the following results for a 7 points finite differences 056 * for abscissae in the [-10, 10] range: 057 * <ul> 058 * <li>step size = 0.25, second order derivative error about 9.97e-10</li> 059 * <li>step size = 0.25, fourth order derivative error about 5.43e-8</li> 060 * <li>step size = 1.0e-6, second order derivative error about 148</li> 061 * <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li> 062 * </ul> 063 * <p> 064 * This example shows that the small step size is really bad, even simply 065 * for second order derivative!</p> 066 * 067 * @since 3.1 068 */ 069public class FiniteDifferencesDifferentiator 070 implements UnivariateFunctionDifferentiator, 071 UnivariateVectorFunctionDifferentiator, 072 UnivariateMatrixFunctionDifferentiator { 073 /** Number of points to use. */ 074 private final int nbPoints; 075 076 /** Step size. */ 077 private final double stepSize; 078 079 /** Half sample span. */ 080 private final double halfSampleSpan; 081 082 /** Lower bound for independent variable. */ 083 private final double tMin; 084 085 /** Upper bound for independent variable. */ 086 private final double tMax; 087 088 /** 089 * Build a differentiator with number of points and step size when independent variable is unbounded. 090 * <p> 091 * Beware that wrong settings for the finite differences differentiator 092 * can lead to highly unstable and inaccurate results, especially for 093 * high derivation orders. Using very small step sizes is often a 094 * <em>bad</em> idea. 095 * </p> 096 * @param nbPoints number of points to use 097 * @param stepSize step size (gap between each point) 098 * @exception NotPositiveException if {@code stepsize <= 0} (note that 099 * {@link NotPositiveException} extends {@link NumberIsTooSmallException}) 100 * @exception NumberIsTooSmallException {@code nbPoint <= 1} 101 */ 102 public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize) 103 throws NotPositiveException, NumberIsTooSmallException { 104 this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY); 105 } 106 107 /** 108 * Build a differentiator with number of points and step size when independent variable is bounded. 109 * <p> 110 * When the independent variable is bounded (tLower < t < tUpper), the sampling 111 * points used for differentiation will be adapted to ensure the constraint holds 112 * even near the boundaries. This means the sample will not be centered anymore in 113 * these cases. At an extreme case, computing derivatives exactly at the lower bound 114 * will lead the sample to be entirely on the right side of the derivation point. 115 * </p> 116 * <p> 117 * Note that the boundaries are considered to be excluded for function evaluation. 118 * </p> 119 * <p> 120 * Beware that wrong settings for the finite differences differentiator 121 * can lead to highly unstable and inaccurate results, especially for 122 * high derivation orders. Using very small step sizes is often a 123 * <em>bad</em> idea. 124 * </p> 125 * @param nbPoints number of points to use 126 * @param stepSize step size (gap between each point) 127 * @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY} 128 * if there are no lower bounds) 129 * @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY} 130 * if there are no upper bounds) 131 * @exception NotPositiveException if {@code stepsize <= 0} (note that 132 * {@link NotPositiveException} extends {@link NumberIsTooSmallException}) 133 * @exception NumberIsTooSmallException {@code nbPoint <= 1} 134 * @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower} 135 */ 136 public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize, 137 final double tLower, final double tUpper) 138 throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 139 140 if (nbPoints <= 1) { 141 throw new NumberIsTooSmallException(stepSize, 1, false); 142 } 143 this.nbPoints = nbPoints; 144 145 if (stepSize <= 0) { 146 throw new NotPositiveException(stepSize); 147 } 148 this.stepSize = stepSize; 149 150 halfSampleSpan = 0.5 * stepSize * (nbPoints - 1); 151 if (2 * halfSampleSpan >= tUpper - tLower) { 152 throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false); 153 } 154 final double safety = JdkMath.ulp(halfSampleSpan); 155 this.tMin = tLower + halfSampleSpan + safety; 156 this.tMax = tUpper - halfSampleSpan - safety; 157 } 158 159 /** 160 * Get the number of points to use. 161 * @return number of points to use 162 */ 163 public int getNbPoints() { 164 return nbPoints; 165 } 166 167 /** 168 * Get the step size. 169 * @return step size 170 */ 171 public double getStepSize() { 172 return stepSize; 173 } 174 175 /** 176 * Evaluate derivatives from a sample. 177 * <p> 178 * Evaluation is done using divided differences. 179 * </p> 180 * @param t evaluation abscissa value and derivatives 181 * @param t0 first sample point abscissa 182 * @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)} 183 * @return value and derivatives at {@code t} 184 * @exception NumberIsTooLargeException if the requested derivation order 185 * is larger or equal to the number of points 186 */ 187 private DerivativeStructure evaluate(final DerivativeStructure t, final double t0, 188 final double[] y) 189 throws NumberIsTooLargeException { 190 191 // create divided differences diagonal arrays 192 final double[] top = new double[nbPoints]; 193 final double[] bottom = new double[nbPoints]; 194 195 for (int i = 0; i < nbPoints; ++i) { 196 197 // update the bottom diagonal of the divided differences array 198 bottom[i] = y[i]; 199 for (int j = 1; j <= i; ++j) { 200 bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize); 201 } 202 203 // update the top diagonal of the divided differences array 204 top[i] = bottom[0]; 205 } 206 207 // evaluate interpolation polynomial (represented by top diagonal) at t 208 final int order = t.getOrder(); 209 final int parameters = t.getFreeParameters(); 210 final double[] derivatives = t.getAllDerivatives(); 211 final double dt0 = t.getValue() - t0; 212 DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0); 213 DerivativeStructure monomial = null; 214 for (int i = 0; i < nbPoints; ++i) { 215 if (i == 0) { 216 // start with monomial(t) = 1 217 monomial = new DerivativeStructure(parameters, order, 1.0); 218 } else { 219 // monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1)) 220 derivatives[0] = dt0 - (i - 1) * stepSize; 221 final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives); 222 monomial = monomial.multiply(deltaX); 223 } 224 interpolation = interpolation.add(monomial.multiply(top[i])); 225 } 226 227 return interpolation; 228 } 229 230 /** {@inheritDoc} 231 * <p>The returned object cannot compute derivatives to arbitrary orders. The 232 * value function will throw a {@link NumberIsTooLargeException} if the requested 233 * derivation order is larger or equal to the number of points. 234 * </p> 235 */ 236 @Override 237 public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) { 238 return new UnivariateDifferentiableFunction() { 239 240 /** {@inheritDoc} */ 241 @Override 242 public double value(final double x) throws MathIllegalArgumentException { 243 return function.value(x); 244 } 245 246 /** {@inheritDoc} */ 247 @Override 248 public DerivativeStructure value(final DerivativeStructure t) 249 throws MathIllegalArgumentException { 250 251 // check we can achieve the requested derivation order with the sample 252 if (t.getOrder() >= nbPoints) { 253 throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); 254 } 255 256 // compute sample position, trying to be centered if possible 257 final double t0 = JdkMath.max(JdkMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; 258 259 // compute sample points 260 final double[] y = new double[nbPoints]; 261 for (int i = 0; i < nbPoints; ++i) { 262 y[i] = function.value(t0 + i * stepSize); 263 } 264 265 // evaluate derivatives 266 return evaluate(t, t0, y); 267 } 268 }; 269 } 270 271 /** {@inheritDoc} 272 * <p>The returned object cannot compute derivatives to arbitrary orders. The 273 * value function will throw a {@link NumberIsTooLargeException} if the requested 274 * derivation order is larger or equal to the number of points. 275 * </p> 276 */ 277 @Override 278 public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) { 279 return new UnivariateDifferentiableVectorFunction() { 280 281 /** {@inheritDoc} */ 282 @Override 283 public double[]value(final double x) throws MathIllegalArgumentException { 284 return function.value(x); 285 } 286 287 /** {@inheritDoc} */ 288 @Override 289 public DerivativeStructure[] value(final DerivativeStructure t) 290 throws MathIllegalArgumentException { 291 292 // check we can achieve the requested derivation order with the sample 293 if (t.getOrder() >= nbPoints) { 294 throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); 295 } 296 297 // compute sample position, trying to be centered if possible 298 final double t0 = JdkMath.max(JdkMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; 299 300 // compute sample points 301 double[][] y = null; 302 for (int i = 0; i < nbPoints; ++i) { 303 final double[] v = function.value(t0 + i * stepSize); 304 if (i == 0) { 305 y = new double[v.length][nbPoints]; 306 } 307 for (int j = 0; j < v.length; ++j) { 308 y[j][i] = v[j]; 309 } 310 } 311 312 // evaluate derivatives 313 final DerivativeStructure[] value = new DerivativeStructure[y.length]; 314 for (int j = 0; j < value.length; ++j) { 315 value[j] = evaluate(t, t0, y[j]); 316 } 317 318 return value; 319 } 320 }; 321 } 322 323 /** {@inheritDoc} 324 * <p>The returned object cannot compute derivatives to arbitrary orders. The 325 * value function will throw a {@link NumberIsTooLargeException} if the requested 326 * derivation order is larger or equal to the number of points. 327 * </p> 328 */ 329 @Override 330 public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) { 331 return new UnivariateDifferentiableMatrixFunction() { 332 333 /** {@inheritDoc} */ 334 @Override 335 public double[][] value(final double x) throws MathIllegalArgumentException { 336 return function.value(x); 337 } 338 339 /** {@inheritDoc} */ 340 @Override 341 public DerivativeStructure[][] value(final DerivativeStructure t) 342 throws MathIllegalArgumentException { 343 344 // check we can achieve the requested derivation order with the sample 345 if (t.getOrder() >= nbPoints) { 346 throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false); 347 } 348 349 // compute sample position, trying to be centered if possible 350 final double t0 = JdkMath.max(JdkMath.min(t.getValue(), tMax), tMin) - halfSampleSpan; 351 352 // compute sample points 353 double[][][] y = null; 354 for (int i = 0; i < nbPoints; ++i) { 355 final double[][] v = function.value(t0 + i * stepSize); 356 if (i == 0) { 357 y = new double[v.length][v[0].length][nbPoints]; 358 } 359 for (int j = 0; j < v.length; ++j) { 360 for (int k = 0; k < v[j].length; ++k) { 361 y[j][k][i] = v[j][k]; 362 } 363 } 364 } 365 366 // evaluate derivatives 367 final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length]; 368 for (int j = 0; j < value.length; ++j) { 369 for (int k = 0; k < y[j].length; ++k) { 370 value[j][k] = evaluate(t, t0, y[j][k]); 371 } 372 } 373 374 return value; 375 } 376 }; 377 } 378}