Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Preliminary checkin of SoC code. 

Contributed by: Xiaogang Zhang


git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@278634 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Joerg Pietschmann 2005-09-04 22:00:27 +00:00
parent a888d20dc6
commit d147d1ebdd
8 changed files with 1133 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

137
src/test/org/apache/commons/math/analysis/DividedDifferenceInterpolatorTest.java Executable file
View File

@ -0,0 +1,137 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Divided Difference interpolator.
* <p>
* The error of polynomial interpolation is
* f(z) - p(z) = f^(n)(zeta) * (z-x[0])(z-x[1])...(z-x[n-1]) / n!
* where f^(n) is the n-th derivative of the approximated function and
* zeta is some point in the interval determined by x[] and z.
* <p>
* Since zeta is unknown, f^(n)(zeta) cannot be calculated. But we can bound
* it and use the absolute value upper bound for estimates. For reference,
* see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 2.
*
* @version $Revision$ $Date$
*/
public final class DividedDifferenceInterpolatorTest extends TestCase {
/**
* Test of interpolator for the sine function.
* <p>
* |sin^(n)(zeta)| <= 1.0, zeta in [0, 2*PI]
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
double x[], y[], z, expected, result, tolerance;
// 6 interpolating points on interval [0, 2*PI]
int n = 6;
double min = 0.0, max = 2 * Math.PI;
x = new double[n];
y = new double[n];
for (int i = 0; i < n; i++) {
x[i] = min + i * (max - min) / n;
y[i] = f.value(x[i]);
}
double derivativebound = 1.0;
UnivariateRealFunction p = interpolator.interpolate(x, y);
z = Math.PI / 4; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = Math.PI * 1.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
}
/**
* Test of interpolator for the exponential function.
* <p>
* |expm1^(n)(zeta)| <= e, zeta in [-1, 1]
*/
public void testExpm1Function() throws MathException {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
double x[], y[], z, expected, result, tolerance;
// 5 interpolating points on interval [-1, 1]
int n = 5;
double min = -1.0, max = 1.0;
x = new double[n];
y = new double[n];
for (int i = 0; i < n; i++) {
x[i] = min + i * (max - min) / n;
y[i] = f.value(x[i]);
}
double derivativebound = Math.E;
UnivariateRealFunction p = interpolator.interpolate(x, y);
z = 0.0; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = 0.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = -0.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the interpolator.
*/
public void testParameters() throws Exception {
UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
try {
// bad abscissas array
double x[] = { 1.0, 2.0, 2.0, 4.0 };
double y[] = { 0.0, 4.0, 4.0, 2.5 };
UnivariateRealFunction p = interpolator.interpolate(x, y);
p.value(0.0);
fail("Expecting MathException - bad abscissas array");
} catch (MathException ex) {
// expected
}
}
/**
* Returns the partial error term (z-x[0])(z-x[1])...(z-x[n-1])/n!
*/
protected double partialerror(double x[], double z) throws
IllegalArgumentException {
if (x.length < 1) {
throw new IllegalArgumentException
("Interpolation array cannot be empty.");
}
double out = 1;
for (int i = 0; i < x.length; i++) {
out *= (z - x[i]) / (i + 1);
}
return out;
}
}

39
src/test/org/apache/commons/math/analysis/Expm1Function.java Executable file
View File

@ -0,0 +1,39 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.FunctionEvaluationException;
/**
* Auxillary class for testing purposes.
*
* @version $Revision$ $Date$
*/
public class Expm1Function implements DifferentiableUnivariateRealFunction {
public double value(double x) throws FunctionEvaluationException {
// Math.expm1() is available in jdk 1.5 but not in jdk 1.4.2.
return Math.exp(x) - 1.0;
}
public UnivariateRealFunction derivative() {
return new UnivariateRealFunction() {
public double value(double x) throws FunctionEvaluationException {
return Math.exp(x);
}
};
}
}

178
src/test/org/apache/commons/math/analysis/LaguerreSolverTest.java Executable file
View File

@ -0,0 +1,178 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import org.apache.commons.math.complex.Complex;
import junit.framework.TestCase;
/**
* Testcase for Laguerre solver.
* <p>
* Laguerre's method is very efficient in solving polynomials. Test runs
* show that for a default absolute accuracy of 1E-6, it generally takes
* less than 5 iterations to find one root, provided solveAll() is not
* invoked, and 15 to 20 iterations to find all roots for quintic function.
*
* @version $Revision$ $Date$
*/
public final class LaguerreSolverTest extends TestCase {
/**
* Test of solver for the linear function.
*/
public void testLinearFunction() throws MathException {
double min, max, expected, result, tolerance;
// p(x) = 4x - 1
double coefficients[] = { -1.0, 4.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver(f);
min = 0.0; max = 1.0; expected = 0.25;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quadratic function.
*/
public void testQuadraticFunction() throws MathException {
double min, max, expected, result, tolerance;
// p(x) = 2x^2 + 5x - 3 = (x+3)(2x-1)
double coefficients[] = { -3.0, 5.0, 2.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver(f);
min = 0.0; max = 2.0; expected = 0.5;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -4.0; max = -1.0; expected = -3.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
double min, max, expected, result, tolerance;
// p(x) = x^5 - x^4 - 12x^3 + x^2 - x - 12 = (x+1)(x+3)(x-4)(x^2-x+1)
double coefficients[] = { -12.0, -1.0, 1.0, -12.0, -1.0, 1.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver(f);
min = -2.0; max = 2.0; expected = -1.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -5.0; max = -2.5; expected = -3.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = 3.0; max = 6.0; expected = 4.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function using solveAll().
*/
public void testQuinticFunction2() throws MathException {
double initial = 0.0, tolerance;
Complex expected, result[];
// p(x) = x^5 + 4x^3 + x^2 + 4 = (x+1)(x^2-x+1)(x^2+4)
double coefficients[] = { 4.0, 0.0, 1.0, 4.0, 0.0, 1.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
LaguerreSolver solver = new LaguerreSolver(f);
result = solver.solveAll(coefficients, initial);
// The order of roots returned by solveAll() depends on
// initial value, solveAll() does no sorting.
expected = new Complex(0.0, -2.0);
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected.abs() * solver.getRelativeAccuracy()));
assertEquals(0.0, (expected.subtract(result[0])).abs(), tolerance);
expected = new Complex(0.0, 2.0);
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected.abs() * solver.getRelativeAccuracy()));
assertEquals(0.0, (expected.subtract(result[1])).abs(), tolerance);
expected = new Complex(0.5, 0.5 * Math.sqrt(3.0));
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected.abs() * solver.getRelativeAccuracy()));
assertEquals(0.0, (expected.subtract(result[2])).abs(), tolerance);
expected = new Complex(-1.0, 0.0);
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected.abs() * solver.getRelativeAccuracy()));
assertEquals(0.0, (expected.subtract(result[3])).abs(), tolerance);
expected = new Complex(0.5, -0.5 * Math.sqrt(3.0));
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected.abs() * solver.getRelativeAccuracy()));
assertEquals(0.0, (expected.subtract(result[4])).abs(), tolerance);
}
/**
* Test of parameters for the solver.
*/
public void testParameters() throws Exception {
double coefficients[] = { -3.0, 5.0, 2.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver(f);
try {
// bad interval
solver.solve(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// no bracketing
solver.solve(2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad function
UnivariateRealFunction f2 = new SinFunction();
UnivariateRealSolver solver2 = new LaguerreSolver(f2);
fail("Expecting IllegalArgumentException - bad function");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

214
src/test/org/apache/commons/math/analysis/MullerSolverTest.java Executable file
View File

@ -0,0 +1,214 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Muller solver.
* <p>
* Muller's method converges almost quadratically near roots, but it can
* be very slow in regions far away from zeros. Test runs show that for
* reasonably good initial values, for a default absolute accuracy of 1E-6,
* it generally takes 5 to 10 iterations for the solver to converge.
* <p>
* Tests for the exponential function illustrate the situations where
* Muller solver performs poorly.
*
* @version $Revision$ $Date$
*/
public final class MullerSolverTest extends TestCase {
/**
* Test of solver for the sine function.
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = Math.PI;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the sine function using solve2().
*/
public void testSinFunction2() throws MathException {
UnivariateRealFunction f = new SinFunction();
MullerSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = Math.PI;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function using solve2().
*/
public void testQuinticFunction2() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
MullerSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function.
* <p>
* It takes 10 to 15 iterations for the last two tests to converge.
* In fact, if not for the bisection alternative, the solver would
* exceed the default maximal iteration of 100.
*/
public void testExpm1Function() throws MathException {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function using solve2().
* <p>
* It takes 25 to 50 iterations for the last two tests to converge.
*/
public void testExpm1Function2() throws MathException {
UnivariateRealFunction f = new Expm1Function();
MullerSolver solver = new MullerSolver(f);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the solver.
*/
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver(f);
try {
// bad interval
solver.solve(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// no bracketing
solver.solve(2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

137
src/test/org/apache/commons/math/analysis/NevilleInterpolatorTest.java Executable file
View File

@ -0,0 +1,137 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Neville interpolator.
* <p>
* The error of polynomial interpolation is
* f(z) - p(z) = f^(n)(zeta) * (z-x[0])(z-x[1])...(z-x[n-1]) / n!
* where f^(n) is the n-th derivative of the approximated function and
* zeta is some point in the interval determined by x[] and z.
* <p>
* Since zeta is unknown, f^(n)(zeta) cannot be calculated. But we can bound
* it and use the absolute value upper bound for estimates. For reference,
* see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 2.
*
* @version $Revision$ $Date$
*/
public final class NevilleInterpolatorTest extends TestCase {
/**
* Test of interpolator for the sine function.
* <p>
* |sin^(n)(zeta)| <= 1.0, zeta in [0, 2*PI]
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealInterpolator interpolator = new NevilleInterpolator();
double x[], y[], z, expected, result, tolerance;
// 6 interpolating points on interval [0, 2*PI]
int n = 6;
double min = 0.0, max = 2 * Math.PI;
x = new double[n];
y = new double[n];
for (int i = 0; i < n; i++) {
x[i] = min + i * (max - min) / n;
y[i] = f.value(x[i]);
}
double derivativebound = 1.0;
UnivariateRealFunction p = interpolator.interpolate(x, y);
z = Math.PI / 4; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = Math.PI * 1.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
}
/**
* Test of interpolator for the exponential function.
* <p>
* |expm1^(n)(zeta)| <= e, zeta in [-1, 1]
*/
public void testExpm1Function() throws MathException {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealInterpolator interpolator = new NevilleInterpolator();
double x[], y[], z, expected, result, tolerance;
// 5 interpolating points on interval [-1, 1]
int n = 5;
double min = -1.0, max = 1.0;
x = new double[n];
y = new double[n];
for (int i = 0; i < n; i++) {
x[i] = min + i * (max - min) / n;
y[i] = f.value(x[i]);
}
double derivativebound = Math.E;
UnivariateRealFunction p = interpolator.interpolate(x, y);
z = 0.0; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = 0.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
z = -0.5; expected = f.value(z); result = p.value(z);
tolerance = Math.abs(derivativebound * partialerror(x, z));
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the interpolator.
*/
public void testParameters() throws Exception {
UnivariateRealInterpolator interpolator = new NevilleInterpolator();
try {
// bad abscissas array
double x[] = { 1.0, 2.0, 2.0, 4.0 };
double y[] = { 0.0, 4.0, 4.0, 2.5 };
UnivariateRealFunction p = interpolator.interpolate(x, y);
p.value(0.0);
fail("Expecting MathException - bad abscissas array");
} catch (MathException ex) {
// expected
}
}
/**
* Returns the partial error term (z-x[0])(z-x[1])...(z-x[n-1])/n!
*/
protected double partialerror(double x[], double z) throws
IllegalArgumentException {
if (x.length < 1) {
throw new IllegalArgumentException
("Interpolation array cannot be empty.");
}
double out = 1;
for (int i = 0; i < x.length; i++) {
out *= (z - x[i]) / (i + 1);
}
return out;
}
}

149
src/test/org/apache/commons/math/analysis/PolynomialFunctionLagrangeFormTest.java Executable file
View File

@ -0,0 +1,149 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Lagrange form of polynomial function.
* <p>
* We use n+1 points to interpolate a polynomial of degree n. This should
* give us the exact same polynomial as result. Thus we can use a very
* small tolerance to account only for round-off errors.
*
* @version $Revision$ $Date$
*/
public final class PolynomialFunctionLagrangeFormTest extends TestCase {
/**
* Test of polynomial for the linear function.
*/
public void testLinearFunction() throws MathException {
PolynomialFunctionLagrangeForm p;
double c[], z, expected, result, tolerance = 1E-12;
// p(x) = 1.5x - 4
double x[] = { 0.0, 3.0 };
double y[] = { -4.0, 0.5 };
p = new PolynomialFunctionLagrangeForm(x, y);
z = 2.0; expected = -1.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 4.5; expected = 2.75; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 6.0; expected = 5.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(1, p.degree());
c = p.getCoefficients();
assertEquals(2, c.length);
assertEquals(-4.0, c[0], tolerance);
assertEquals(1.5, c[1], tolerance);
}
/**
* Test of polynomial for the quadratic function.
*/
public void testQuadraticFunction() throws MathException {
PolynomialFunctionLagrangeForm p;
double c[], z, expected, result, tolerance = 1E-12;
// p(x) = 2x^2 + 5x - 3 = (2x - 1)(x + 3)
double x[] = { 0.0, -1.0, 0.5 };
double y[] = { -3.0, -6.0, 0.0 };
p = new PolynomialFunctionLagrangeForm(x, y);
z = 1.0; expected = 4.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 2.5; expected = 22.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = -2.0; expected = -5.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(2, p.degree());
c = p.getCoefficients();
assertEquals(3, c.length);
assertEquals(-3.0, c[0], tolerance);
assertEquals(5.0, c[1], tolerance);
assertEquals(2.0, c[2], tolerance);
}
/**
* Test of polynomial for the quintic function.
*/
public void testQuinticFunction() throws MathException {
PolynomialFunctionLagrangeForm p;
double c[], z, expected, result, tolerance = 1E-12;
// p(x) = x^5 - x^4 - 7x^3 + x^2 + 6x = x(x^2 - 1)(x + 2)(x - 3)
double x[] = { 1.0, -1.0, 2.0, 3.0, -3.0, 0.5 };
double y[] = { 0.0, 0.0, -24.0, 0.0, -144.0, 2.34375 };
p = new PolynomialFunctionLagrangeForm(x, y);
z = 0.0; expected = 0.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = -2.0; expected = 0.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 4.0; expected = 360.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(5, p.degree());
c = p.getCoefficients();
assertEquals(6, c.length);
assertEquals(0.0, c[0], tolerance);
assertEquals(6.0, c[1], tolerance);
assertEquals(1.0, c[2], tolerance);
assertEquals(-7.0, c[3], tolerance);
assertEquals(-1.0, c[4], tolerance);
assertEquals(1.0, c[5], tolerance);
}
/**
* Test of parameters for the polynomial.
*/
public void testParameters() throws Exception {
PolynomialFunctionLagrangeForm p;
try {
// bad input array length
double x[] = { 1.0 };
double y[] = { 2.0 };
p = new PolynomialFunctionLagrangeForm(x, y);
fail("Expecting IllegalArgumentException - bad input array length");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// mismatch input arrays
double x[] = { 1.0, 2.0, 3.0, 4.0 };
double y[] = { 0.0, -4.0, -24.0 };
p = new PolynomialFunctionLagrangeForm(x, y);
fail("Expecting IllegalArgumentException - mismatch input arrays");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

148
src/test/org/apache/commons/math/analysis/PolynomialFunctionNewtonFormTest.java Executable file
View File

@ -0,0 +1,148 @@
/*
* Copyright 2003-2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Newton form of polynomial function.
* <p>
* The small tolerance number is used only to account for round-off errors.
*
* @version $Revision$ $Date$
*/
public final class PolynomialFunctionNewtonFormTest extends TestCase {
/**
* Test of polynomial for the linear function.
*/
public void testLinearFunction() throws MathException {
PolynomialFunctionNewtonForm p;
double coefficients[], z, expected, result, tolerance = 1E-12;
// p(x) = 1.5x - 4 = 2 + 1.5(x-4)
double a[] = { 2.0, 1.5 };
double c[] = { 4.0 };
p = new PolynomialFunctionNewtonForm(a, c);
z = 2.0; expected = -1.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 4.5; expected = 2.75; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 6.0; expected = 5.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(1, p.degree());
coefficients = p.getCoefficients();
assertEquals(2, coefficients.length);
assertEquals(-4.0, coefficients[0], tolerance);
assertEquals(1.5, coefficients[1], tolerance);
}
/**
* Test of polynomial for the quadratic function.
*/
public void testQuadraticFunction() throws MathException {
PolynomialFunctionNewtonForm p;
double coefficients[], z, expected, result, tolerance = 1E-12;
// p(x) = 2x^2 + 5x - 3 = 4 + 3(x-1) + 2(x-1)(x+2)
double a[] = { 4.0, 3.0, 2.0 };
double c[] = { 1.0, -2.0 };
p = new PolynomialFunctionNewtonForm(a, c);
z = 1.0; expected = 4.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 2.5; expected = 22.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = -2.0; expected = -5.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(2, p.degree());
coefficients = p.getCoefficients();
assertEquals(3, coefficients.length);
assertEquals(-3.0, coefficients[0], tolerance);
assertEquals(5.0, coefficients[1], tolerance);
assertEquals(2.0, coefficients[2], tolerance);
}
/**
* Test of polynomial for the quintic function.
*/
public void testQuinticFunction() throws MathException {
PolynomialFunctionNewtonForm p;
double coefficients[], z, expected, result, tolerance = 1E-12;
// p(x) = x^5 - x^4 - 7x^3 + x^2 + 6x
// = 6x - 6x^2 -6x^2(x-1) + x^2(x-1)(x+1) + x^2(x-1)(x+1)(x-2)
double a[] = { 0.0, 6.0, -6.0, -6.0, 1.0, 1.0 };
double c[] = { 0.0, 0.0, 1.0, -1.0, 2.0 };
p = new PolynomialFunctionNewtonForm(a, c);
z = 0.0; expected = 0.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = -2.0; expected = 0.0; result = p.value(z);
assertEquals(expected, result, tolerance);
z = 4.0; expected = 360.0; result = p.value(z);
assertEquals(expected, result, tolerance);
assertEquals(5, p.degree());
coefficients = p.getCoefficients();
assertEquals(6, coefficients.length);
assertEquals(0.0, coefficients[0], tolerance);
assertEquals(6.0, coefficients[1], tolerance);
assertEquals(1.0, coefficients[2], tolerance);
assertEquals(-7.0, coefficients[3], tolerance);
assertEquals(-1.0, coefficients[4], tolerance);
assertEquals(1.0, coefficients[5], tolerance);
}
/**
* Test of parameters for the polynomial.
*/
public void testParameters() throws Exception {
PolynomialFunctionNewtonForm p;
try {
// bad input array length
double a[] = { 1.0 };
double c[] = { 2.0 };
p = new PolynomialFunctionNewtonForm(a, c);
fail("Expecting IllegalArgumentException - bad input array length");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// mismatch input arrays
double a[] = { 1.0, 2.0, 3.0, 4.0 };
double c[] = { 4.0, 3.0, 2.0, 1.0 };
p = new PolynomialFunctionNewtonForm(a, c);
fail("Expecting IllegalArgumentException - mismatch input arrays");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

131
src/test/org/apache/commons/math/analysis/RiddersSolverTest.java Executable file
View File

@ -0,0 +1,131 @@
/*
* Copyright 2005 The Apache Software Foundation.
*
* Licensed 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.math.analysis;
import org.apache.commons.math.MathException;
import junit.framework.TestCase;
/**
* Testcase for Ridders solver.
* <p>
* Ridders' method converges superlinearly, more specific, its rate of
* convergence is sqrt(2). Test runs show that for a default absolute
* accuracy of 1E-6, it generally takes less than 5 iterations for close
* initial bracket and 5 to 10 iterations for distant initial bracket
* to converge.
*
* @version $Revision$ $Date$
*/
public final class RiddersSolverTest extends TestCase {
/**
* Test of solver for the sine function.
*/
public void testSinFunction() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new RiddersSolver(f);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = Math.PI;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealSolver solver = new RiddersSolver(f);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function.
*/
public void testExpm1Function() throws MathException {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealSolver solver = new RiddersSolver(f);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = Math.max(solver.getAbsoluteAccuracy(),
Math.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the solver.
*/
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new RiddersSolver(f);
try {
// bad interval
solver.solve(1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// no bracketing
solver.solve(2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
}
}