mirror of
https://github.com/apache/commons-math.git
synced 2025-03-03 15:09:09 +00:00
Added Hermite interpolator.
This class implements the Hermite polynomial interpolation method, which can match function derivatives in addition to function value at sampling points. The implementation is done for vector-valued functions. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1351257 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
5f00c9b777
commit
e8a6708cfd
@ -54,7 +54,13 @@ All rights reserved
|
||||
The LocalizedFormatsTest class in the unit tests is an adapted version of
|
||||
the OrekitMessagesTest class from the orekit library distributed under the
|
||||
terms of the Apache 2 licence. Original source copyright:
|
||||
Copyright 2010 CS Communication & Systèmes
|
||||
Copyright 2010 CS Systèmes d'Information
|
||||
===============================================================================
|
||||
|
||||
The HermiteInterpolator class and its corresponding test have been imported from
|
||||
the orekit library distributed under the terms of the Apache 2 licence. Original
|
||||
source copyright:
|
||||
Copyright 2010-2012 CS Systèmes d'Information
|
||||
===============================================================================
|
||||
|
||||
The complete text of licenses and disclaimers associated with the the original
|
||||
|
||||
@ -52,6 +52,11 @@ If the output is not quite correct, check for invisible trailing spaces!
|
||||
<body>
|
||||
<release version="3.1" date="TBD" description="
|
||||
">
|
||||
<action dev="luc" type="add">
|
||||
A new HermiteInterpolator class allow interpolation of vector-valued
|
||||
functions using both values and derivatives of the function at sample
|
||||
points.
|
||||
</action>
|
||||
<action dev="erans" type="fix" issue="MATH-804">
|
||||
Parameterized "CurveFitter" class (package "o.a.c.m.optimization.fitting")
|
||||
with the type of the fitting function. Updated subclasses "PolynomialFitter",
|
||||
|
||||
@ -0,0 +1,254 @@
|
||||
/*
|
||||
* 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.interpolation;
|
||||
|
||||
import java.util.ArrayList;
|
||||
import java.util.List;
|
||||
|
||||
import org.apache.commons.math3.analysis.DifferentiableUnivariateVectorFunction;
|
||||
import org.apache.commons.math3.analysis.UnivariateVectorFunction;
|
||||
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
|
||||
import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
||||
import org.apache.commons.math3.exception.MathIllegalStateException;
|
||||
import org.apache.commons.math3.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math3.util.ArithmeticUtils;
|
||||
|
||||
/** Polynomial interpolator using both sample values and sample derivatives.
|
||||
* <p>
|
||||
* The interpolation polynomials match all sample points, including both values
|
||||
* and provided derivatives. There is one polynomial for each component of
|
||||
* the values vector. All polynomial have the same degree. The degree of the
|
||||
* polynomials depends on the number of points and number of derivatives at each
|
||||
* point. For example the interpolation polynomials for n sample points without
|
||||
* any derivatives all have degree n-1. The interpolation polynomials for n
|
||||
* sample points with the two extreme points having value and first derivative
|
||||
* and the remaining points having value only all have degree n+1. The
|
||||
* interpolation polynomial for n sample points with value, first and second
|
||||
* derivative for all points all have degree 3n-1.
|
||||
* </p>
|
||||
* <p>
|
||||
* This class has been imported from the Orekit space flight dynamics library
|
||||
* also distributed under the terms of the Apache License V2. Original copyright
|
||||
* is: Copyright 2002-2012 CS Systèmes d'Information.
|
||||
* </p>
|
||||
* @version $Id$
|
||||
* @since 3.1
|
||||
*/
|
||||
public class HermiteInterpolator implements DifferentiableUnivariateVectorFunction {
|
||||
|
||||
/** Sample abscissae. */
|
||||
private final List<Double> abscissae;
|
||||
|
||||
/** Top diagonal of the divided differences array. */
|
||||
private final List<double[]> topDiagonal;
|
||||
|
||||
/** Bottom diagonal of the divided differences array. */
|
||||
private final List<double[]> bottomDiagonal;
|
||||
|
||||
/** Create an empty interpolator.
|
||||
*/
|
||||
public HermiteInterpolator() {
|
||||
this.abscissae = new ArrayList<Double>();
|
||||
this.topDiagonal = new ArrayList<double[]>();
|
||||
this.bottomDiagonal = new ArrayList<double[]>();
|
||||
}
|
||||
|
||||
/** Add a sample point.
|
||||
* <p>
|
||||
* This method must be called once for each sample point. It is allowed to
|
||||
* mix some calls with values only with calls with values and first
|
||||
* derivatives.
|
||||
* </p>
|
||||
* <p>
|
||||
* The point abscissae for all calls <em>must</em> be different.
|
||||
* </p>
|
||||
* @param x abscissa of the sample point
|
||||
* @param value value and derivatives of the sample point
|
||||
* (if only one row is passed, it is the value, if two rows are
|
||||
* passed the first one is the value and the second the derivative
|
||||
* and so on)
|
||||
* @exception MathIllegalArgumentException if the abscissa is equals to a previously
|
||||
* added sample point
|
||||
*/
|
||||
public void addSamplePoint(final double x, final double[] ... value)
|
||||
throws MathIllegalArgumentException {
|
||||
|
||||
for (int i = 0; i < value.length; ++i) {
|
||||
|
||||
final double[] y = value[i].clone();
|
||||
if (i > 1) {
|
||||
double inv = 1.0 / ArithmeticUtils.factorial(i);
|
||||
for (int j = 0; j < y.length; ++j) {
|
||||
y[j] *= inv;
|
||||
}
|
||||
}
|
||||
|
||||
// update the bottom diagonal of the divided differences array
|
||||
final int n = abscissae.size();
|
||||
bottomDiagonal.add(n - i, y);
|
||||
double[] bottom0 = y;
|
||||
for (int j = i; j < n; ++j) {
|
||||
final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
|
||||
final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
|
||||
if (Double.isInfinite(inv)) {
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO,
|
||||
x);
|
||||
}
|
||||
for (int k = 0; k < y.length; ++k) {
|
||||
bottom1[k] = inv * (bottom0[k] - bottom1[k]);
|
||||
}
|
||||
bottom0 = bottom1;
|
||||
}
|
||||
|
||||
// update the top diagonal of the divided differences array
|
||||
topDiagonal.add(bottom0.clone());
|
||||
|
||||
// update the abscissae array
|
||||
abscissae.add(x);
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Compute the interpolation polynomials.
|
||||
* @return interpolation polynomials array
|
||||
* @exception MathIllegalStateException if sample is empty
|
||||
*/
|
||||
public PolynomialFunction[] getPolynomials()
|
||||
throws MathIllegalStateException {
|
||||
|
||||
// safety check
|
||||
checkInterpolation();
|
||||
|
||||
// iteration initialization
|
||||
final PolynomialFunction zero = polynomial(0);
|
||||
PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
|
||||
for (int i = 0; i < polynomials.length; ++i) {
|
||||
polynomials[i] = zero;
|
||||
}
|
||||
PolynomialFunction coeff = polynomial(1);
|
||||
|
||||
// build the polynomials by iterating on the top diagonal of the divided differences array
|
||||
for (int i = 0; i < topDiagonal.size(); ++i) {
|
||||
double[] tdi = topDiagonal.get(i);
|
||||
for (int k = 0; k < polynomials.length; ++k) {
|
||||
polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
|
||||
}
|
||||
coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
|
||||
}
|
||||
|
||||
return polynomials;
|
||||
|
||||
}
|
||||
|
||||
/** Interpolate value at a specified abscissa.
|
||||
* <p>
|
||||
* Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
|
||||
* value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
|
||||
* except it does not build the intermediate polynomials, so this method is faster and
|
||||
* numerically more stable.
|
||||
* </p>
|
||||
* @param x interpolation abscissa
|
||||
* @return interpolated value
|
||||
* @exception MathIllegalStateException if sample is empty
|
||||
*/
|
||||
public double[] value(double x)
|
||||
throws MathIllegalStateException {
|
||||
|
||||
// safety check
|
||||
checkInterpolation();
|
||||
|
||||
final double[] value = new double[topDiagonal.get(0).length];
|
||||
double valueCoeff = 1;
|
||||
for (int i = 0; i < topDiagonal.size(); ++i) {
|
||||
double[] dividedDifference = topDiagonal.get(i);
|
||||
for (int k = 0; k < value.length; ++k) {
|
||||
value[k] += dividedDifference[k] * valueCoeff;
|
||||
}
|
||||
final double deltaX = x - abscissae.get(i);
|
||||
valueCoeff *= deltaX;
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
}
|
||||
|
||||
/** Interpolate first derivative at a specified abscissa.
|
||||
* <p>
|
||||
* Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
|
||||
* value} methods of the derivatives of all polynomials returned by {@link
|
||||
* #getPolynomials() getPolynomials}, except it builds neither the intermediate
|
||||
* polynomials nor their derivatives, so this method is faster and numerically more stable.
|
||||
* </p>
|
||||
* @param x interpolation abscissa
|
||||
* @return interpolated derivative
|
||||
* @exception MathIllegalStateException if sample is empty
|
||||
*/
|
||||
public double[] derivative(double x)
|
||||
throws MathIllegalStateException {
|
||||
|
||||
// safety check
|
||||
checkInterpolation();
|
||||
|
||||
final double[] derivative = new double[topDiagonal.get(0).length];
|
||||
double valueCoeff = 1;
|
||||
double derivativeCoeff = 0;
|
||||
for (int i = 0; i < topDiagonal.size(); ++i) {
|
||||
double[] dividedDifference = topDiagonal.get(i);
|
||||
for (int k = 0; k < derivative.length; ++k) {
|
||||
derivative[k] += dividedDifference[k] * derivativeCoeff;
|
||||
}
|
||||
final double deltaX = x - abscissae.get(i);
|
||||
derivativeCoeff = valueCoeff + derivativeCoeff * deltaX;
|
||||
valueCoeff *= deltaX;
|
||||
}
|
||||
|
||||
return derivative;
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc}} */
|
||||
public UnivariateVectorFunction derivative() {
|
||||
return new UnivariateVectorFunction() {
|
||||
|
||||
/** {@inheritDoc}} */
|
||||
public double[] value(double x) {
|
||||
return derivative(x);
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
/** Check interpolation can be performed.
|
||||
* @exception MathIllegalStateException if interpolation cannot be performed
|
||||
* because sample is empty
|
||||
*/
|
||||
private void checkInterpolation() throws MathIllegalStateException {
|
||||
if (abscissae.isEmpty()) {
|
||||
throw new MathIllegalStateException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
|
||||
}
|
||||
}
|
||||
|
||||
/** Create a polynomial from its coefficients.
|
||||
* @param c polynomials coefficients
|
||||
* @return polynomial
|
||||
*/
|
||||
private PolynomialFunction polynomial(double ... c) {
|
||||
return new PolynomialFunction(c);
|
||||
}
|
||||
|
||||
}
|
||||
@ -88,9 +88,10 @@ public enum LocalizedFormats implements Localizable {
|
||||
DIMENSIONS_MISMATCH("dimensions mismatch"), /* keep */
|
||||
DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN("Discrete cumulative probability function returned NaN for argument {0}"),
|
||||
DISTRIBUTION_NOT_LOADED("distribution not loaded"),
|
||||
DUPLICATED_ABSCISSA("Abscissa {0} is duplicated at both indices {1} and {2}"),
|
||||
DUPLICATED_ABSCISSA_DIVISION_BY_ZERO("duplicated abscissa {0} causes division by zero"),
|
||||
ELITISM_RATE("elitism rate ({0})"),
|
||||
EMPTY_CLUSTER_IN_K_MEANS("empty cluster in k-means"),
|
||||
EMPTY_INTERPOLATION_SAMPLE("sample for interpolation is empty"),
|
||||
EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY("empty polynomials coefficients array"), /* keep */
|
||||
EMPTY_SELECTED_COLUMN_INDEX_ARRAY("empty selected column index array"),
|
||||
EMPTY_SELECTED_ROW_INDEX_ARRAY("empty selected row index array"),
|
||||
@ -112,7 +113,6 @@ public enum LocalizedFormats implements Localizable {
|
||||
GCD_OVERFLOW_32_BITS("overflow: gcd({0}, {1}) is 2^31"),
|
||||
GCD_OVERFLOW_64_BITS("overflow: gcd({0}, {1}) is 2^63"),
|
||||
HOLE_BETWEEN_MODELS_TIME_RANGES("{0} wide hole between models time ranges"),
|
||||
IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO("identical abscissas x[{0}] == x[{1}] == {2} cause division by zero"),
|
||||
ILL_CONDITIONED_OPERATOR("condition number {1} is too high "),
|
||||
INDEX_LARGER_THAN_MAX("the index specified: {0} is larger than the current maximal index {1}"),
|
||||
INDEX_NOT_POSITIVE("index ({0}) is not positive"),
|
||||
|
||||
@ -60,9 +60,10 @@ DIMENSIONS_MISMATCH_SIMPLE = {0} != {1}
|
||||
DIMENSIONS_MISMATCH = dimensions incoh\u00e9rentes
|
||||
DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN = Discr\u00e8tes fonction de probabilit\u00e9 cumulative retourn\u00e9 NaN \u00e0 l''argument de {0}
|
||||
DISTRIBUTION_NOT_LOADED = aucune distribution n''a \u00e9t\u00e9 charg\u00e9e
|
||||
DUPLICATED_ABSCISSA = Abscisse {0} dupliqu\u00e9e aux indices {1} et {2}
|
||||
DUPLICATED_ABSCISSA_DIVISION_BY_ZERO = la duplication de l''abscisse {0} engendre une division par z\u00e9ro
|
||||
ELITISM_RATE = proportion d''\u00e9litisme ({0})
|
||||
EMPTY_CLUSTER_IN_K_MEANS = groupe vide dans l''algorithme des k-moyennes
|
||||
EMPTY_INTERPOLATION_SAMPLE = \u00e9chantillon d''interpolation vide
|
||||
EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY = tableau de coefficients polyn\u00f4miaux vide
|
||||
EMPTY_SELECTED_COLUMN_INDEX_ARRAY = tableau des indices de colonnes s\u00e9lectionn\u00e9es vide
|
||||
EMPTY_SELECTED_ROW_INDEX_ARRAY = tableau des indices de lignes s\u00e9lectionn\u00e9es vide
|
||||
@ -84,7 +85,6 @@ FUNCTION_NOT_POLYNOMIAL = la fonction n''est pas p\u00f4lynomiale
|
||||
GCD_OVERFLOW_32_BITS = d\u00e9passement de capacit\u00e9 : le PGCD de {0} et {1} vaut 2^31
|
||||
GCD_OVERFLOW_64_BITS = d\u00e9passement de capacit\u00e9 : le PGCD de {0} et {1} vaut 2^63
|
||||
HOLE_BETWEEN_MODELS_TIME_RANGES = trou de longueur {0} entre les domaines temporels des mod\u00e8les
|
||||
IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO = les abscisses identiques x[{0}] == x[{1}] == {2} engendrent une division par z\u00e9ro
|
||||
ILL_CONDITIONED_OPERATOR = le conditionnement {1} est trop \u00e9lev\u00e9
|
||||
INDEX_LARGER_THAN_MAX = l''index sp\u00e9cifi\u00e9 ({0}) d\u00e9passe l''index maximal courant ({1})
|
||||
INDEX_NOT_POSITIVE = l''indice ({0}) n''est pas positif
|
||||
|
||||
@ -436,6 +436,27 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
|
||||
It has been described in William Dudziak's <a
|
||||
href="http://www.dudziak.com/microsphere.pdf">MS thesis</a>.
|
||||
</p>
|
||||
<p>
|
||||
<a href="http://en.wikipedia.org/wiki/Hermite_interpolation">Hermite interpolation</a>
|
||||
is an interpolation method that can use derivatives in addition to function values at sample points. The <a
|
||||
href="../apidocs/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.html">HermiteInterpolator</a>
|
||||
class implements this method for vector-valued functions. The sampling points can have any spacing (there are
|
||||
no requirements for a regular grid) and some points may provide derivatives while others don't provide them
|
||||
(or provide derivatives to a smaller order). Points are added one at a time, as shown in the following example:
|
||||
</p>
|
||||
<source>HermiteInterpolator interpolator = new HermiteInterpolator;
|
||||
// at x = 0, we provide both value and first derivative
|
||||
interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
|
||||
// at x = 1, we provide only function value
|
||||
interpolator.addSamplePoint(1.0, new double[] { 4.0 });
|
||||
// at x = 2, we provide both value and first derivative
|
||||
interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
|
||||
// should print "value at x = 0.5: 2.5625"
|
||||
System.out.println("value at x = 0.5: " + interpolator.value(0.5)[0]);
|
||||
// should print "derivative at x = 0.5: 3.5"
|
||||
System.out.println("derivative at x = 0.5: " + interpolator.derivative(0.5)[0]);
|
||||
// should print "interpolation polynomial: 1 + 2 x + 4 x^2 - 4 x^3 + x^4"
|
||||
System.out.println("interpolation polynomial: " + interpolator.getPolynomials()[0]);</source>
|
||||
<p>
|
||||
A <a href="../apidocs/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.html">
|
||||
BivariateGridInterpolator</a> is used to find a bivariate real-valued
|
||||
|
||||
@ -0,0 +1,246 @@
|
||||
/*
|
||||
* 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.interpolation;
|
||||
|
||||
import java.util.Random;
|
||||
|
||||
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
||||
public class HermiteInterpolatorTest {
|
||||
|
||||
@Test
|
||||
public void testZero() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(0.0, new double[] { 0.0 });
|
||||
for (double x = -10; x < 10; x += 1.0) {
|
||||
Assert.assertEquals(0.0, interpolator.value(x)[0], 1.0e-15);
|
||||
Assert.assertEquals(0.0, interpolator.derivative(x)[0], 1.0e-15);
|
||||
}
|
||||
checkPolynomial(new PolynomialFunction(new double[] { 0.0 }),
|
||||
interpolator.getPolynomials()[0]);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testQuadratic() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(0.0, new double[] { 2.0 });
|
||||
interpolator.addSamplePoint(1.0, new double[] { 0.0 });
|
||||
interpolator.addSamplePoint(2.0, new double[] { 0.0 });
|
||||
for (double x = -10; x < 10; x += 1.0) {
|
||||
Assert.assertEquals((x - 1.0) * (x - 2.0), interpolator.value(x)[0], 1.0e-15);
|
||||
Assert.assertEquals(2 * x - 3.0, interpolator.derivative(x)[0], 1.0e-15);
|
||||
}
|
||||
checkPolynomial(new PolynomialFunction(new double[] { 2.0, -3.0, 1.0 }),
|
||||
interpolator.getPolynomials()[0]);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testMixedDerivatives() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
|
||||
interpolator.addSamplePoint(1.0, new double[] { 4.0 });
|
||||
interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
|
||||
Assert.assertEquals(4, interpolator.getPolynomials()[0].degree());
|
||||
Assert.assertEquals(1.0, interpolator.value(0.0)[0], 1.0e-15);
|
||||
Assert.assertEquals(2.0, interpolator.derivative(0.0)[0], 1.0e-15);
|
||||
Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15);
|
||||
Assert.assertEquals(5.0, interpolator.value(2.0)[0], 1.0e-15);
|
||||
Assert.assertEquals(2.0, interpolator.derivative(2.0)[0], 1.0e-15);
|
||||
checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }),
|
||||
interpolator.getPolynomials()[0]);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testRandomPolynomialsValuesOnly() {
|
||||
|
||||
Random random = new Random(0x42b1e7dbd361a932l);
|
||||
|
||||
for (int i = 0; i < 100; ++i) {
|
||||
|
||||
int maxDegree = 0;
|
||||
PolynomialFunction[] p = new PolynomialFunction[5];
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
int degree = random.nextInt(7);
|
||||
p[k] = randomPolynomial(degree, random);
|
||||
maxDegree = FastMath.max(maxDegree, degree);
|
||||
}
|
||||
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
for (int j = 0; j < 1 + maxDegree; ++j) {
|
||||
double x = 0.1 * j;
|
||||
double[] values = new double[p.length];
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
values[k] = p[k].value(x);
|
||||
}
|
||||
interpolator.addSamplePoint(x, values);
|
||||
}
|
||||
|
||||
for (double x = 0; x < 2; x += 0.1) {
|
||||
double[] values = interpolator.value(x);
|
||||
Assert.assertEquals(p.length, values.length);
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * FastMath.abs(p[k].value(x)));
|
||||
}
|
||||
}
|
||||
|
||||
PolynomialFunction[] result = interpolator.getPolynomials();
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
checkPolynomial(p[k], result[k]);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testRandomPolynomialsFirstDerivative() {
|
||||
|
||||
Random random = new Random(0x570803c982ca5d3bl);
|
||||
|
||||
for (int i = 0; i < 100; ++i) {
|
||||
|
||||
int maxDegree = 0;
|
||||
PolynomialFunction[] p = new PolynomialFunction[5];
|
||||
PolynomialFunction[] pPrime = new PolynomialFunction[5];
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
int degree = random.nextInt(7);
|
||||
p[k] = randomPolynomial(degree, random);
|
||||
pPrime[k] = p[k].polynomialDerivative();
|
||||
maxDegree = FastMath.max(maxDegree, degree);
|
||||
}
|
||||
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
for (int j = 0; j < 1 + maxDegree / 2; ++j) {
|
||||
double x = 0.1 * j;
|
||||
double[] values = new double[p.length];
|
||||
double[] derivatives = new double[p.length];
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
values[k] = p[k].value(x);
|
||||
derivatives[k] = pPrime[k].value(x);
|
||||
}
|
||||
interpolator.addSamplePoint(x, values, derivatives);
|
||||
}
|
||||
|
||||
for (double x = 0; x < 2; x += 0.1) {
|
||||
double[] values = interpolator.value(x);
|
||||
double[] derivatives = interpolator.derivative(x);
|
||||
Assert.assertEquals(p.length, values.length);
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * FastMath.abs(p[k].value(x)));
|
||||
Assert.assertEquals(pPrime[k].value(x), derivatives[k], 4.0e-8 * FastMath.abs(p[k].value(x)));
|
||||
}
|
||||
}
|
||||
|
||||
PolynomialFunction[] result = interpolator.getPolynomials();
|
||||
for (int k = 0; k < p.length; ++k) {
|
||||
checkPolynomial(p[k], result[k]);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testSine() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
for (double x = 0; x < FastMath.PI; x += 0.5) {
|
||||
interpolator.addSamplePoint(x, new double[] { FastMath.sin(x) });
|
||||
}
|
||||
for (double x = 0.1; x <= 2.9; x += 0.01) {
|
||||
Assert.assertEquals(FastMath.sin(x), interpolator.value(x)[0], 3.5e-5);
|
||||
Assert.assertEquals(FastMath.cos(x), interpolator.derivative(x)[0], 1.3e-4);
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testSquareRoot() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
for (double x = 1.0; x < 3.6; x += 0.5) {
|
||||
interpolator.addSamplePoint(x, new double[] { FastMath.sqrt(x) });
|
||||
}
|
||||
for (double x = 1.1; x < 3.5; x += 0.01) {
|
||||
Assert.assertEquals(FastMath.sqrt(x), interpolator.value(x)[0], 1.5e-4);
|
||||
Assert.assertEquals(0.5 / FastMath.sqrt(x), interpolator.derivative(x)[0], 8.5e-4);
|
||||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testWikipedia() {
|
||||
// this test corresponds to the example from Wikipedia page:
|
||||
// http://en.wikipedia.org/wiki/Hermite_interpolation
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 });
|
||||
interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] { 0 }, new double[] { 0 });
|
||||
interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] { 8 }, new double[] { 56 });
|
||||
for (double x = -1.0; x <= 1.0; x += 0.125) {
|
||||
double x2 = x * x;
|
||||
double x4 = x2 * x2;
|
||||
double x8 = x4 * x4;
|
||||
Assert.assertEquals(x8 + 1, interpolator.value(x)[0], 1.0e-15);
|
||||
Assert.assertEquals(8 * x4 * x2 * x, interpolator.derivative(x)[0], 1.0e-15);
|
||||
}
|
||||
checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }),
|
||||
interpolator.getPolynomials()[0]);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testOnePointParabola() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(0, new double[] { 1 }, new double[] { 1 }, new double[] { 2 });
|
||||
for (double x = -1.0; x <= 1.0; x += 0.125) {
|
||||
Assert.assertEquals(1 + x * (1 + x), interpolator.value(x)[0], 1.0e-15);
|
||||
Assert.assertEquals(1 + 2 * x, interpolator.derivative(x)[0], 1.0e-15);
|
||||
}
|
||||
checkPolynomial(new PolynomialFunction(new double[] { 1, 1, 1 }),
|
||||
interpolator.getPolynomials()[0]);
|
||||
}
|
||||
|
||||
private PolynomialFunction randomPolynomial(int degree, Random random) {
|
||||
double[] coeff = new double[ 1 + degree];
|
||||
for (int j = 0; j < degree; ++j) {
|
||||
coeff[j] = random.nextDouble();
|
||||
}
|
||||
return new PolynomialFunction(coeff);
|
||||
}
|
||||
|
||||
@Test(expected=IllegalStateException.class)
|
||||
public void testEmptySample() {
|
||||
new HermiteInterpolator().value(0.0);
|
||||
}
|
||||
|
||||
@Test(expected=IllegalArgumentException.class)
|
||||
public void testDuplicatedAbscissa() {
|
||||
HermiteInterpolator interpolator = new HermiteInterpolator();
|
||||
interpolator.addSamplePoint(1.0, new double[] { 0.0 });
|
||||
interpolator.addSamplePoint(1.0, new double[] { 1.0 });
|
||||
}
|
||||
|
||||
private void checkPolynomial(PolynomialFunction expected, PolynomialFunction result) {
|
||||
Assert.assertTrue(result.degree() >= expected.degree());
|
||||
double[] cE = expected.getCoefficients();
|
||||
double[] cR = result.getCoefficients();
|
||||
for (int i = 0; i < cE.length; ++i) {
|
||||
Assert.assertEquals(cE[i], cR[i], 1.0e-8 * FastMath.abs(cE[i]));
|
||||
}
|
||||
for (int i = cE.length; i < cR.length; ++i) {
|
||||
Assert.assertEquals(0.0, cR[i], 1.0e-9);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user