Added Hermite interpolator for ExtendFieldElement instances.
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1449657 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
49d94ad78c
commit
684a87be70
|
|
@ -55,6 +55,9 @@ This is a minor release: It combines bug fixes and new features.
|
||||||
Changes to existing features were made in a backwards-compatible
|
Changes to existing features were made in a backwards-compatible
|
||||||
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
|
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
|
||||||
">
|
">
|
||||||
|
<action dev="luc" type="add" >
|
||||||
|
Added Hermite interpolator for ExtendFieldElement instances.
|
||||||
|
</action>
|
||||||
<action dev="luc" type="add" >
|
<action dev="luc" type="add" >
|
||||||
Added ExtendFieldElement interface to represent anything that is
|
Added ExtendFieldElement interface to represent anything that is
|
||||||
real number like, implemented by both Decimal64, Dfp and DerivativeStructure.
|
real number like, implemented by both Decimal64, Dfp and DerivativeStructure.
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,152 @@
|
||||||
|
/*
|
||||||
|
* 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.FieldElement;
|
||||||
|
import org.apache.commons.math3.exception.MathArithmeticException;
|
||||||
|
import org.apache.commons.math3.exception.NoDataException;
|
||||||
|
import org.apache.commons.math3.exception.ZeroException;
|
||||||
|
import org.apache.commons.math3.exception.util.LocalizedFormats;
|
||||||
|
import org.apache.commons.math3.util.MathArrays;
|
||||||
|
|
||||||
|
/** 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 polynomials 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>
|
||||||
|
*
|
||||||
|
* @version $Id$
|
||||||
|
* @since 3.2
|
||||||
|
*/
|
||||||
|
public class FieldHermiteInterpolator<T extends FieldElement<T>> {
|
||||||
|
|
||||||
|
/** Sample abscissae. */
|
||||||
|
private final List<T> abscissae;
|
||||||
|
|
||||||
|
/** Top diagonal of the divided differences array. */
|
||||||
|
private final List<T[]> topDiagonal;
|
||||||
|
|
||||||
|
/** Bottom diagonal of the divided differences array. */
|
||||||
|
private final List<T[]> bottomDiagonal;
|
||||||
|
|
||||||
|
/** Create an empty interpolator.
|
||||||
|
*/
|
||||||
|
public FieldHermiteInterpolator() {
|
||||||
|
this.abscissae = new ArrayList<T>();
|
||||||
|
this.topDiagonal = new ArrayList<T[]>();
|
||||||
|
this.bottomDiagonal = new ArrayList<T[]>();
|
||||||
|
}
|
||||||
|
|
||||||
|
/** 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 ZeroException if the abscissa difference between added point
|
||||||
|
* and a previous point is zero (i.e. the two points are at same abscissa)
|
||||||
|
* @exception MathArithmeticException if the number of derivatives is larger
|
||||||
|
* than 20, which prevents computation of a factorial
|
||||||
|
*/
|
||||||
|
public void addSamplePoint(final T x, final T[] ... value)
|
||||||
|
throws ZeroException, MathArithmeticException {
|
||||||
|
|
||||||
|
T factorial = x.getField().getOne();
|
||||||
|
for (int i = 0; i < value.length; ++i) {
|
||||||
|
|
||||||
|
final T[] y = value[i].clone();
|
||||||
|
if (i > 1) {
|
||||||
|
factorial = factorial.multiply(i);
|
||||||
|
final T inv = factorial.reciprocal();
|
||||||
|
for (int j = 0; j < y.length; ++j) {
|
||||||
|
y[j] = y[j].multiply(inv);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// update the bottom diagonal of the divided differences array
|
||||||
|
final int n = abscissae.size();
|
||||||
|
bottomDiagonal.add(n - i, y);
|
||||||
|
T[] bottom0 = y;
|
||||||
|
for (int j = i; j < n; ++j) {
|
||||||
|
final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
|
||||||
|
if (x.equals(abscissae.get(n - (j + 1)))) {
|
||||||
|
throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
|
||||||
|
}
|
||||||
|
final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal();
|
||||||
|
for (int k = 0; k < y.length; ++k) {
|
||||||
|
bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
|
||||||
|
}
|
||||||
|
bottom0 = bottom1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// update the top diagonal of the divided differences array
|
||||||
|
topDiagonal.add(bottom0.clone());
|
||||||
|
|
||||||
|
// update the abscissae array
|
||||||
|
abscissae.add(x);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/** Interpolate value at a specified abscissa.
|
||||||
|
* @param x interpolation abscissa
|
||||||
|
* @return interpolated value
|
||||||
|
* @exception NoDataException if sample is empty
|
||||||
|
*/
|
||||||
|
public T[] value(T x) throws NoDataException {
|
||||||
|
|
||||||
|
// safety check
|
||||||
|
if (abscissae.isEmpty()) {
|
||||||
|
throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
|
||||||
|
}
|
||||||
|
|
||||||
|
final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
|
||||||
|
T valueCoeff = x.getField().getOne();
|
||||||
|
for (int i = 0; i < topDiagonal.size(); ++i) {
|
||||||
|
T[] dividedDifference = topDiagonal.get(i);
|
||||||
|
for (int k = 0; k < value.length; ++k) {
|
||||||
|
value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff));
|
||||||
|
}
|
||||||
|
final T deltaX = x.subtract(abscissae.get(i));
|
||||||
|
valueCoeff = valueCoeff.multiply(deltaX);
|
||||||
|
}
|
||||||
|
|
||||||
|
return value;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,245 @@
|
||||||
|
/*
|
||||||
|
* 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.dfp.Dfp;
|
||||||
|
import org.apache.commons.math3.dfp.DfpField;
|
||||||
|
import org.apache.commons.math3.exception.NoDataException;
|
||||||
|
import org.apache.commons.math3.fraction.BigFraction;
|
||||||
|
import org.apache.commons.math3.util.FastMath;
|
||||||
|
import org.junit.Assert;
|
||||||
|
import org.junit.Test;
|
||||||
|
|
||||||
|
public class FieldHermiteInterpolatorTest {
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testZero() {
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] { new BigFraction(0) });
|
||||||
|
for (int x = -10; x < 10; x++) {
|
||||||
|
BigFraction y = interpolator.value(new BigFraction(x))[0];
|
||||||
|
Assert.assertEquals(BigFraction.ZERO, y);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testQuadratic() {
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] { new BigFraction(2) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] { new BigFraction(0) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction(2), new BigFraction[] { new BigFraction(0) });
|
||||||
|
for (double x = -10; x < 10; x += 1.0) {
|
||||||
|
BigFraction y = interpolator.value(new BigFraction(x))[0];
|
||||||
|
Assert.assertEquals((x - 1) * (x - 2), y.doubleValue(), 1.0e-15);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testMixedDerivatives() {
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] { new BigFraction(1) }, new BigFraction[] { new BigFraction(2) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] { new BigFraction(4) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction(2), new BigFraction[] { new BigFraction(5) }, new BigFraction[] { new BigFraction(2) });
|
||||||
|
Assert.assertEquals(new BigFraction(1), interpolator.value(new BigFraction(0))[0]);
|
||||||
|
Assert.assertEquals(new BigFraction(4), interpolator.value(new BigFraction(1))[0]);
|
||||||
|
Assert.assertEquals(new BigFraction(5), interpolator.value(new BigFraction(2))[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);
|
||||||
|
}
|
||||||
|
|
||||||
|
DfpField field = new DfpField(30);
|
||||||
|
Dfp step = field.getOne().divide(field.newDfp(10));
|
||||||
|
FieldHermiteInterpolator<Dfp> interpolator = new FieldHermiteInterpolator<Dfp>();
|
||||||
|
for (int j = 0; j < 1 + maxDegree; ++j) {
|
||||||
|
Dfp x = field.newDfp(j).multiply(step);
|
||||||
|
Dfp[] values = new Dfp[p.length];
|
||||||
|
for (int k = 0; k < p.length; ++k) {
|
||||||
|
values[k] = field.newDfp(p[k].value(x.getReal()));
|
||||||
|
}
|
||||||
|
interpolator.addSamplePoint(x, values);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int j = 0; j < 20; ++j) {
|
||||||
|
Dfp x = field.newDfp(j).multiply(step);
|
||||||
|
Dfp[] values = interpolator.value(x);
|
||||||
|
Assert.assertEquals(p.length, values.length);
|
||||||
|
for (int k = 0; k < p.length; ++k) {
|
||||||
|
Assert.assertEquals(p[k].value(x.getReal()),
|
||||||
|
values[k].getReal(),
|
||||||
|
1.0e-8 * FastMath.abs(p[k].value(x.getReal())));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
@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);
|
||||||
|
}
|
||||||
|
|
||||||
|
DfpField field = new DfpField(30);
|
||||||
|
Dfp step = field.getOne().divide(field.newDfp(10));
|
||||||
|
FieldHermiteInterpolator<Dfp> interpolator = new FieldHermiteInterpolator<Dfp>();
|
||||||
|
for (int j = 0; j < 1 + maxDegree / 2; ++j) {
|
||||||
|
Dfp x = field.newDfp(j).multiply(step);
|
||||||
|
Dfp[] values = new Dfp[p.length];
|
||||||
|
Dfp[] derivatives = new Dfp[p.length];
|
||||||
|
for (int k = 0; k < p.length; ++k) {
|
||||||
|
values[k] = field.newDfp(p[k].value(x.getReal()));
|
||||||
|
derivatives[k] = field.newDfp(pPrime[k].value(x.getReal()));
|
||||||
|
}
|
||||||
|
interpolator.addSamplePoint(x, values, derivatives);
|
||||||
|
}
|
||||||
|
|
||||||
|
Dfp h = step.divide(field.newDfp(100000));
|
||||||
|
for (int j = 0; j < 20; ++j) {
|
||||||
|
Dfp x = field.newDfp(j).multiply(step);
|
||||||
|
Dfp[] y = interpolator.value(x);
|
||||||
|
Dfp[] yP = interpolator.value(x.add(h));
|
||||||
|
Dfp[] yM = interpolator.value(x.subtract(h));
|
||||||
|
Assert.assertEquals(p.length, y.length);
|
||||||
|
for (int k = 0; k < p.length; ++k) {
|
||||||
|
Assert.assertEquals(p[k].value(x.getReal()),
|
||||||
|
y[k].getReal(),
|
||||||
|
1.0e-8 * FastMath.abs(p[k].value(x.getReal())));
|
||||||
|
Assert.assertEquals(pPrime[k].value(x.getReal()),
|
||||||
|
yP[k].subtract(yM[k]).divide(h.multiply(2)).getReal(),
|
||||||
|
4.0e-8 * FastMath.abs(p[k].value(x.getReal())));
|
||||||
|
}
|
||||||
|
System.out.println();
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testSine() {
|
||||||
|
DfpField field = new DfpField(30);
|
||||||
|
FieldHermiteInterpolator<Dfp> interpolator = new FieldHermiteInterpolator<Dfp>();
|
||||||
|
for (Dfp x = field.getZero(); x.getReal() < FastMath.PI; x = x.add(0.5)) {
|
||||||
|
interpolator.addSamplePoint(x, new Dfp[] { x.sin() });
|
||||||
|
}
|
||||||
|
for (Dfp x = field.newDfp(0.1); x.getReal() < 2.9; x = x.add(0.01)) {
|
||||||
|
Dfp y = interpolator.value(x)[0];
|
||||||
|
Assert.assertEquals( x.sin().getReal(), y.getReal(), 3.5e-5);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testSquareRoot() {
|
||||||
|
DfpField field = new DfpField(30);
|
||||||
|
FieldHermiteInterpolator<Dfp> interpolator = new FieldHermiteInterpolator<Dfp>();
|
||||||
|
for (Dfp x = field.getOne(); x.getReal() < 3.6; x = x.add(0.5)) {
|
||||||
|
interpolator.addSamplePoint(x, new Dfp[] { x.sqrt() });
|
||||||
|
}
|
||||||
|
for (Dfp x = field.newDfp(1.1); x.getReal() < 3.5; x = x.add(0.01)) {
|
||||||
|
Dfp y = interpolator.value(x)[0];
|
||||||
|
Assert.assertEquals(x.sqrt().getReal(), y.getReal(), 1.5e-4);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testWikipedia() {
|
||||||
|
// this test corresponds to the example from Wikipedia page:
|
||||||
|
// http://en.wikipedia.org/wiki/Hermite_interpolation
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(-1),
|
||||||
|
new BigFraction[] { new BigFraction( 2) },
|
||||||
|
new BigFraction[] { new BigFraction(-8) },
|
||||||
|
new BigFraction[] { new BigFraction(56) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction( 0),
|
||||||
|
new BigFraction[] { new BigFraction( 1) },
|
||||||
|
new BigFraction[] { new BigFraction( 0) },
|
||||||
|
new BigFraction[] { new BigFraction( 0) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction( 1),
|
||||||
|
new BigFraction[] { new BigFraction( 2) },
|
||||||
|
new BigFraction[] { new BigFraction( 8) },
|
||||||
|
new BigFraction[] { new BigFraction(56) });
|
||||||
|
for (BigFraction x = new BigFraction(-1); x.doubleValue() <= 1.0; x = x.add(new BigFraction(1, 8))) {
|
||||||
|
BigFraction y = interpolator.value(x)[0];
|
||||||
|
BigFraction x2 = x.multiply(x);
|
||||||
|
BigFraction x4 = x2.multiply(x2);
|
||||||
|
BigFraction x8 = x4.multiply(x4);
|
||||||
|
Assert.assertEquals(x8.add(new BigFraction(1)), y);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
public void testOnePointParabola() {
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(0),
|
||||||
|
new BigFraction[] { new BigFraction(1) },
|
||||||
|
new BigFraction[] { new BigFraction(1) },
|
||||||
|
new BigFraction[] { new BigFraction(2) });
|
||||||
|
for (BigFraction x = new BigFraction(-1); x.doubleValue() <= 1.0; x = x.add(new BigFraction(1, 8))) {
|
||||||
|
BigFraction y = interpolator.value(x)[0];
|
||||||
|
Assert.assertEquals(BigFraction.ONE.add(x.multiply(BigFraction.ONE.add(x))), y);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
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=NoDataException.class)
|
||||||
|
public void testEmptySample() {
|
||||||
|
new FieldHermiteInterpolator<BigFraction>().value(BigFraction.ZERO);
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test(expected=IllegalArgumentException.class)
|
||||||
|
public void testDuplicatedAbscissa() {
|
||||||
|
FieldHermiteInterpolator<BigFraction> interpolator = new FieldHermiteInterpolator<BigFraction>();
|
||||||
|
interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] { new BigFraction(0) });
|
||||||
|
interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] { new BigFraction(1) });
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
Loading…
Reference in New Issue