FiniteDifferencesDifferentiator.java

  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.  */

  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;

  24. import java.io.Serializable;

  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;

  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 {

  78.     /** Serializable UID. */
  79.     private static final long serialVersionUID = 20120917L;

  81.     /** Number of points to use. */
  82.     private final int nbPoints;

  84.     /** Step size. */
  85.     private final double stepSize;

  87.     /** Half sample span. */
  88.     private final double halfSampleSpan;

  90.     /** Lower bound for independent variable. */
  91.     private final double tMin;

  93.     /** Upper bound for independent variable. */
  94.     private final double tMax;

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

  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 {

  148.         if (nbPoints <= 1) {
  149.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  150.                                                    stepSize, 1);
  151.         }
  152.         this.nbPoints = nbPoints;

  154.         if (stepSize <= 0) {
  155.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
  156.                                                    stepSize, 0);
  157.         }
  158.         this.stepSize = stepSize;

  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;

  169.     }

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

  179.     /**
  180.      * Get the step size.
  181.      * @return step size
  182.      */
  183.     public double getStepSize() {
  184.         return stepSize;
  185.     }

  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 {

  203.         // create divided differences diagonal arrays
  204.         final double[] top    = new double[nbPoints];
  205.         final double[] bottom = new double[nbPoints];

  207.         for (int i = 0; i < nbPoints; ++i) {

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

  215.             // update the top diagonal of the divided differences array
  216.             top[i] = bottom[0];

  218.         }

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

  235.         return interpolation;

  237.     }

  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() {

  249.             /** {@inheritDoc} */
  250.             @Override
  251.             public double value(final double x) throws MathIllegalArgumentException {
  252.                 return function.value(x);
  253.             }

  255.             /** {@inheritDoc} */
  256.             @Override
  257.             public <T extends Derivative<T>> T value(T t)
  258.                 throws MathIllegalArgumentException {

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

  266.                 // compute sample position, trying to be centered if possible
  267.                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;

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

  275.                 // evaluate derivatives
  276.                 return evaluate(t, t0, y);

  278.             }

  280.         };
  281.     }

  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() {

  293.             /** {@inheritDoc} */
  294.             @Override
  295.             public double[]value(final double x) throws MathIllegalArgumentException {
  296.                 return function.value(x);
  297.             }

  299.             /** {@inheritDoc} */
  300.             @Override
  301.             public <T extends Derivative<T>> T[] value(T t)
  302.                 throws MathIllegalArgumentException {

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

  310.                 // compute sample position, trying to be centered if possible
  311.                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;

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

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

  331.                 return value;

  333.             }

  335.         };
  336.     }

  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() {

  348.             /** {@inheritDoc} */
  349.             @Override
  350.             public double[][]  value(final double x) throws MathIllegalArgumentException {
  351.                 return function.value(x);
  352.             }

  354.             /** {@inheritDoc} */
  355.             @Override
  356.             public <T extends Derivative<T>> T[][] value(T t)
  357.                 throws MathIllegalArgumentException {

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

  365.                 // compute sample position, trying to be centered if possible
  366.                 final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;

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

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

  390.                 return value;

  392.             }

  394.         };
  395.     }

  397. }
Created with JaCoCo 0.8.12.202403310830