Added a FiniteDifferencesDifferentiator class.
This class implements both UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator and UnivariateMatrixFunctionDifferentiator. It is not limited in the derivation order, but checks the order is consistent with the number of points. Typically, attempting to extract 5th derivative from a 3 points scheme will not work (a NumberIsTooLargeException will be thrown, referencing the derivation order). git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1386745 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
3876bd9bc4
commit
830cd9eeca
|
|
@ -0,0 +1,306 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math3.analysis.differentiation;
|
||||
|
||||
import java.io.Serializable;
|
||||
|
||||
import org.apache.commons.math3.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateVectorFunction;
|
||||
import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
||||
import org.apache.commons.math3.exception.NotPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
|
||||
/** Univariate functions differentiator using finite differences.
|
||||
* <p>
|
||||
* This class creates some wrapper objetcs around regular
|
||||
* {@link UnivariateFunction univariate functions} (or {@link
|
||||
* UnivariateVectorFunction univariate vector functions} or {@link
|
||||
* UnivariateMatrixFunction univariate matrix functions}). These
|
||||
* wrapper objects compute derivatives in addition to function
|
||||
* value.
|
||||
* </p>
|
||||
* <p>
|
||||
* The wrapper objects work by calling the underlying function on
|
||||
* a sampling grid around the current point and perform polynomial
|
||||
* interpolation. A finite differences scheme with n points is
|
||||
* theoretically able to compute derivatives up to order n-1, but
|
||||
* it is generally better to have a slight margin. The step size must
|
||||
* also be small enough in order for the polynomial approximation to
|
||||
* be good in the current point neighborhood, but it should not be too
|
||||
* small because numerical instability appears quickly (there are several
|
||||
* differences of close points). Choosing the number of points and
|
||||
* the step size is highly problem dependent.
|
||||
* </p>
|
||||
* <p>
|
||||
* As an example of good and bad settings, lets consider the quintic
|
||||
* polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}.
|
||||
* Since it is a polynomial, finite differences with at least 6 points
|
||||
* should theoretically recover the exact same polynomial and hence
|
||||
* compute accurate derivatives for any order. However, due to numerical
|
||||
* errors, we get the following results for a 7 points finite differences
|
||||
* for abscissae in the [-10, 10] range:
|
||||
* <ul>
|
||||
* <li>step size = 0.25, second order derivative error about 9.97e-10</li>
|
||||
* <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
|
||||
* <li>step size = 1.0e-6, second order derivative error about 56.25</li>
|
||||
* <li>step size = 1.0e-6, fourth order derivative error about 2.47e+14</li>
|
||||
* </ul>
|
||||
* This example shows that the small step size is really bad, even simply
|
||||
* for second order derivative!
|
||||
* </p>
|
||||
* @version $Id$
|
||||
* @since 3.1
|
||||
*/
|
||||
public class FiniteDifferencesDifferentiator
|
||||
implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
|
||||
UnivariateMatrixFunctionDifferentiator, Serializable {
|
||||
|
||||
/** Serializable UID. */
|
||||
private static final long serialVersionUID = 20120917L;
|
||||
|
||||
/** Number of points to use. */
|
||||
private final int nbPoints;
|
||||
|
||||
/** Step size. */
|
||||
private double stepSize;
|
||||
|
||||
/**
|
||||
* Build a differentiator with number of points and step size.
|
||||
* <p>
|
||||
* Beware that wrong settings for the finite differences differentiator
|
||||
* can lead to highly unstable and inaccurate results, especially for
|
||||
* high derivation orders. Using very small step sizes is often a
|
||||
* <em>bad</em> idea.
|
||||
* </p>
|
||||
* @param nbPoints number of points to use
|
||||
* @param stepSize step size (gap between each point)
|
||||
* @exception NotPositiveException if {@code stepsize <= 0} (note that
|
||||
* {@link NotPositiveException} extends {@link NumberIsTooSmallException})
|
||||
* @exception NumberIsTooSmallException {@code nbPoint <= 1}
|
||||
*/
|
||||
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
|
||||
throws NotPositiveException, NumberIsTooSmallException {
|
||||
|
||||
if (nbPoints <= 1) {
|
||||
throw new NumberIsTooSmallException(stepSize, 1, false);
|
||||
}
|
||||
this.nbPoints = nbPoints;
|
||||
|
||||
if (stepSize <= 0) {
|
||||
throw new NotPositiveException(stepSize);
|
||||
}
|
||||
this.stepSize = stepSize;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the number of points to use.
|
||||
* @return number of points to use
|
||||
*/
|
||||
public int getNbPoints() {
|
||||
return nbPoints;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the step size.
|
||||
* @return step size
|
||||
*/
|
||||
public double getStepSize() {
|
||||
return stepSize;
|
||||
}
|
||||
|
||||
/**
|
||||
* Evaluate derivatives from a centered sample.
|
||||
* @param t central value and derivatives
|
||||
* @param y function values at {@code t + stepSize * (i - 0.5 * (nbPoints - 1))}
|
||||
* @return value and derivatives at {@code t}
|
||||
* @exception NumberIsTooLargeException if the requested derivation order
|
||||
* is larger or equal to the number of points
|
||||
*/
|
||||
private DerivativeStructure evaluate(final DerivativeStructure t, final double[] y)
|
||||
throws NumberIsTooLargeException {
|
||||
|
||||
// check we can achieve the requested derivation order with the sample
|
||||
final int order = t.getOrder();
|
||||
if (order >= nbPoints) {
|
||||
throw new NumberIsTooLargeException(order, nbPoints, false);
|
||||
}
|
||||
|
||||
// create divided differences diagonal arrays
|
||||
final double[] top = new double[nbPoints];
|
||||
final double[] bottom = new double[nbPoints];
|
||||
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
|
||||
// update the bottom diagonal of the divided differences array
|
||||
bottom[i] = y[i];
|
||||
for (int j = 1; j <= i; ++j) {
|
||||
bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize);
|
||||
}
|
||||
|
||||
// update the top diagonal of the divided differences array
|
||||
top[i] = bottom[0];
|
||||
|
||||
}
|
||||
|
||||
// evaluate interpolation polynomial (represented by top diagonal) at t
|
||||
final int parameters = t.getFreeParameters();
|
||||
final double[] derivatives = t.getAllDerivatives();
|
||||
DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0);
|
||||
DerivativeStructure monomial = new DerivativeStructure(parameters, order, 1.0);
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
interpolation = interpolation.add(monomial.multiply(top[i]));
|
||||
derivatives[0] = stepSize * (0.5 * (nbPoints - 1) - i);
|
||||
final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
|
||||
monomial = monomial.multiply(deltaX);
|
||||
}
|
||||
|
||||
return interpolation;
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc}
|
||||
* <p>The returned object cannot compute derivatives to arbitrary orders. The
|
||||
* value function will throw a {@link NumberIsTooLargeException} if the requested
|
||||
* derivation order is larger or equal to the number of points.
|
||||
* </p>
|
||||
*/
|
||||
public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) {
|
||||
return new UnivariateDifferentiableFunction() {
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double value(final double x) throws MathIllegalArgumentException {
|
||||
return function.value(x);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public DerivativeStructure value(final DerivativeStructure t)
|
||||
throws MathIllegalArgumentException {
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
final double[] y = new double[nbPoints];
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
y[i] = function.value(xi);
|
||||
}
|
||||
|
||||
// evaluate derivatives
|
||||
return evaluate(t, y);
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
/** {@inheritDoc}
|
||||
* <p>The returned object cannot compute derivatives to arbitrary orders. The
|
||||
* value function will throw a {@link NumberIsTooLargeException} if the requested
|
||||
* derivation order is larger or equal to the number of points.
|
||||
* </p>
|
||||
*/
|
||||
public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) {
|
||||
return new UnivariateDifferentiableVectorFunction() {
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[]value(final double x) throws MathIllegalArgumentException {
|
||||
return function.value(x);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public DerivativeStructure[] value(final DerivativeStructure t)
|
||||
throws MathIllegalArgumentException {
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
double[][] y = null;
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
final double[] v = function.value(xi);
|
||||
if (i == 0) {
|
||||
y = new double[v.length][nbPoints];
|
||||
}
|
||||
for (int j = 0; j < v.length; ++j) {
|
||||
y[j][i] = v[j];
|
||||
}
|
||||
}
|
||||
|
||||
// evaluate derivatives
|
||||
final DerivativeStructure[] value = new DerivativeStructure[y.length];
|
||||
for (int j = 0; j < value.length; ++j) {
|
||||
value[j] = evaluate(t, y[j]);
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
/** {@inheritDoc}
|
||||
* <p>The returned object cannot compute derivatives to arbitrary orders. The
|
||||
* value function will throw a {@link NumberIsTooLargeException} if the requested
|
||||
* derivation order is larger or equal to the number of points.
|
||||
* </p>
|
||||
*/
|
||||
public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) {
|
||||
return new UnivariateDifferentiableMatrixFunction() {
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[][] value(final double x) throws MathIllegalArgumentException {
|
||||
return function.value(x);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public DerivativeStructure[][] value(final DerivativeStructure t)
|
||||
throws MathIllegalArgumentException {
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
double[][][] y = null;
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
final double[][] v = function.value(xi);
|
||||
if (i == 0) {
|
||||
y = new double[v.length][v[0].length][nbPoints];
|
||||
}
|
||||
for (int j = 0; j < v.length; ++j) {
|
||||
for (int k = 0; k < v[j].length; ++k) {
|
||||
y[j][k][i] = v[j][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// evaluate derivatives
|
||||
final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length];
|
||||
for (int j = 0; j < value.length; ++j) {
|
||||
for (int k = 0; k < y[j].length; ++k) {
|
||||
value[j][k] = evaluate(t, y[j][k]);
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,264 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math3.analysis.differentiation;
|
||||
|
||||
import org.apache.commons.math3.TestUtils;
|
||||
import org.apache.commons.math3.analysis.QuinticFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateVectorFunction;
|
||||
import org.apache.commons.math3.analysis.function.Gaussian;
|
||||
import org.apache.commons.math3.analysis.function.Sin;
|
||||
import org.apache.commons.math3.exception.NotPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
||||
/**
|
||||
* Test for class {@link FiniteDifferencesDifferentiator}.
|
||||
*/
|
||||
public class FiniteDifferencesDifferentiatorTest {
|
||||
|
||||
@Test(expected=NumberIsTooSmallException.class)
|
||||
public void testWrongNumberOfPoints() {
|
||||
new FiniteDifferencesDifferentiator(1, 1.0);
|
||||
}
|
||||
|
||||
@Test(expected=NotPositiveException.class)
|
||||
public void testWrongStepSize() {
|
||||
new FiniteDifferencesDifferentiator(3, 0.0);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testSerialization() {
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(3, 1.0e-3);
|
||||
FiniteDifferencesDifferentiator recovered =
|
||||
(FiniteDifferencesDifferentiator) TestUtils.serializeAndRecover(differentiator);
|
||||
Assert.assertEquals(differentiator.getNbPoints(), recovered.getNbPoints());
|
||||
Assert.assertEquals(differentiator.getStepSize(), recovered.getStepSize(), 1.0e-15);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testConstant() {
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(5, 0.01);
|
||||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
return 42.0;
|
||||
}
|
||||
});
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
Assert.assertEquals(42.0, y.getValue(), 1.0e-15);
|
||||
Assert.assertEquals( 0.0, y.getPartialDerivative(1), 1.0e-15);
|
||||
Assert.assertEquals( 0.0, y.getPartialDerivative(2), 1.0e-15);
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testLinear() {
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(5, 0.01);
|
||||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
return 2 - 3 * x;
|
||||
}
|
||||
});
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
Assert.assertEquals(2 - 3 * x, y.getValue(), 1.0e-20);
|
||||
Assert.assertEquals(-3.0, y.getPartialDerivative(1), 4.0e-13);
|
||||
Assert.assertEquals( 0.0, y.getPartialDerivative(2), 5.0e-11);
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testGaussian() {
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(9, 0.02);
|
||||
UnivariateDifferentiableFunction gaussian = new Gaussian(1.0, 2.0);
|
||||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(gaussian);
|
||||
double[] expectedError = new double[] {
|
||||
2.776e-17, 1.742e-15, 2.385e-13, 1.329e-11, 2.668e-9, 8.873e-8
|
||||
};
|
||||
double[] maxError = new double[expectedError.length];
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure dsX = new DerivativeStructure(1, maxError.length - 1, 0, x);
|
||||
DerivativeStructure yRef = gaussian.value(dsX);
|
||||
DerivativeStructure y = f.value(dsX);
|
||||
for (int order = 0; order <= yRef.getOrder(); ++order) {
|
||||
maxError[order] = FastMath.max(maxError[order],
|
||||
FastMath.abs(yRef.getPartialDerivative(order) -
|
||||
y.getPartialDerivative(order)));
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < maxError.length; ++i) {
|
||||
Assert.assertEquals(expectedError[i], maxError[i], 0.01 * expectedError[i]);
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testStepSizeUnstability() {
|
||||
UnivariateDifferentiableFunction quintic = new QuinticFunction();
|
||||
UnivariateDifferentiableFunction goodStep =
|
||||
new FiniteDifferencesDifferentiator(7, 0.25).differentiate(quintic);
|
||||
UnivariateDifferentiableFunction badStep =
|
||||
new FiniteDifferencesDifferentiator(7, 1.0e-6).differentiate(quintic);
|
||||
double[] maxErrorGood = new double[7];
|
||||
double[] maxErrorBad = new double[7];
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure dsX = new DerivativeStructure(1, 6, 0, x);
|
||||
DerivativeStructure yRef = quintic.value(dsX);
|
||||
DerivativeStructure yGood = goodStep.value(dsX);
|
||||
DerivativeStructure yBad = badStep.value(dsX);
|
||||
for (int order = 0; order <= 6; ++order) {
|
||||
maxErrorGood[order] = FastMath.max(maxErrorGood[order],
|
||||
FastMath.abs(yRef.getPartialDerivative(order) -
|
||||
yGood.getPartialDerivative(order)));
|
||||
maxErrorBad[order] = FastMath.max(maxErrorBad[order],
|
||||
FastMath.abs(yRef.getPartialDerivative(order) -
|
||||
yBad.getPartialDerivative(order)));
|
||||
}
|
||||
}
|
||||
|
||||
// the 0.25 step size is good for finite differences in the quintic on this abscissa range for 7 points
|
||||
// the errors are fair
|
||||
final double[] expectedGood = new double[] {
|
||||
7.276e-12, 7.276e-11, 9.968e-10, 3.092e-9, 5.432e-8, 8.196e-8, 1.818e-6
|
||||
};
|
||||
|
||||
// the 1.0e-6 step size is far too small for finite differences in the quintic on this abscissa range for 7 points
|
||||
// the errors are huge!
|
||||
final double[] expectedBad = new double[] {
|
||||
1.792e-22, 6.926e-5, 56.25, 1.783e8, 2.468e14, 3.056e20, 5.857e26
|
||||
};
|
||||
|
||||
for (int i = 0; i < maxErrorGood.length; ++i) {
|
||||
Assert.assertEquals(expectedGood[i], maxErrorGood[i], 0.01 * expectedGood[i]);
|
||||
Assert.assertEquals(expectedBad[i], maxErrorBad[i], 0.01 * expectedBad[i]);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testVectorFunction() {
|
||||
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(7, 0.01);
|
||||
UnivariateDifferentiableVectorFunction f =
|
||||
differentiator.differentiate(new UnivariateVectorFunction() {
|
||||
|
||||
public double[] value(double x) {
|
||||
return new double[] { FastMath.cos(x), FastMath.sin(x) };
|
||||
}
|
||||
|
||||
});
|
||||
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure[] y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
double cos = FastMath.cos(x);
|
||||
double sin = FastMath.sin(x);
|
||||
Assert.assertEquals(cos, y[0].getValue(), 2.0e-16);
|
||||
Assert.assertEquals(sin, y[1].getValue(), 2.0e-16);
|
||||
Assert.assertEquals(-sin, y[0].getPartialDerivative(1), 5.0e-14);
|
||||
Assert.assertEquals( cos, y[1].getPartialDerivative(1), 5.0e-14);
|
||||
Assert.assertEquals(-cos, y[0].getPartialDerivative(2), 6.0e-12);
|
||||
Assert.assertEquals(-sin, y[1].getPartialDerivative(2), 6.0e-12);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testMatrixFunction() {
|
||||
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(7, 0.01);
|
||||
UnivariateDifferentiableMatrixFunction f =
|
||||
differentiator.differentiate(new UnivariateMatrixFunction() {
|
||||
|
||||
public double[][] value(double x) {
|
||||
return new double[][] {
|
||||
{ FastMath.cos(x), FastMath.sin(x) },
|
||||
{ FastMath.cosh(x), FastMath.sinh(x) }
|
||||
};
|
||||
}
|
||||
|
||||
});
|
||||
|
||||
for (double x = -1; x < 1; x += 0.02) {
|
||||
DerivativeStructure[][] y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
double cos = FastMath.cos(x);
|
||||
double sin = FastMath.sin(x);
|
||||
double cosh = FastMath.cosh(x);
|
||||
double sinh = FastMath.sinh(x);
|
||||
Assert.assertEquals(cos, y[0][0].getValue(), 7.0e-18);
|
||||
Assert.assertEquals(sin, y[0][1].getValue(), 7.0e-18);
|
||||
Assert.assertEquals(cosh, y[1][0].getValue(), 3.0e-16);
|
||||
Assert.assertEquals(sinh, y[1][1].getValue(), 3.0e-16);
|
||||
Assert.assertEquals(-sin, y[0][0].getPartialDerivative(1), 2.0e-14);
|
||||
Assert.assertEquals( cos, y[0][1].getPartialDerivative(1), 2.0e-14);
|
||||
Assert.assertEquals( sinh, y[1][0].getPartialDerivative(1), 3.0e-14);
|
||||
Assert.assertEquals( cosh, y[1][1].getPartialDerivative(1), 3.0e-14);
|
||||
Assert.assertEquals(-cos, y[0][0].getPartialDerivative(2), 3.0e-12);
|
||||
Assert.assertEquals(-sin, y[0][1].getPartialDerivative(2), 3.0e-12);
|
||||
Assert.assertEquals( cosh, y[1][0].getPartialDerivative(2), 6.0e-12);
|
||||
Assert.assertEquals( sinh, y[1][1].getPartialDerivative(2), 6.0e-12);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testSeveralFreeParameters() {
|
||||
FiniteDifferencesDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(5, 0.001);
|
||||
UnivariateDifferentiableFunction sine = new Sin();
|
||||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(sine);
|
||||
double[] expectedError = new double[] {
|
||||
1.110e-16, 2.66e-12, 4.803e-9, 5.486e-5
|
||||
};
|
||||
double[] maxError = new double[expectedError.length];
|
||||
for (double x = -2; x < 2; x += 0.1) {
|
||||
for (double y = -2; y < 2; y += 0.1) {
|
||||
DerivativeStructure dsX = new DerivativeStructure(2, maxError.length - 1, 0, x);
|
||||
DerivativeStructure dsY = new DerivativeStructure(2, maxError.length - 1, 1, y);
|
||||
DerivativeStructure dsT = dsX.multiply(3).subtract(dsY.multiply(2));
|
||||
DerivativeStructure sRef = sine.value(dsT);
|
||||
DerivativeStructure s = f.value(dsT);
|
||||
for (int xOrder = 0; xOrder <= sRef.getOrder(); ++xOrder) {
|
||||
for (int yOrder = 0; yOrder <= sRef.getOrder(); ++yOrder) {
|
||||
if (xOrder + yOrder <= sRef.getOrder()) {
|
||||
maxError[xOrder +yOrder] = FastMath.max(maxError[xOrder + yOrder],
|
||||
FastMath.abs(sRef.getPartialDerivative(xOrder, yOrder) -
|
||||
s.getPartialDerivative(xOrder, yOrder)));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < maxError.length; ++i) {
|
||||
Assert.assertEquals(expectedError[i], maxError[i], 0.01 * expectedError[i]);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue