Author: luc
Date: Fri Jul 29 15:44:21 2011
New Revision: 1152276
URL: http://svn.apache.org/viewvc?rev=1152276&view=rev
Log:
Added a Brent-like solver that has higher (user specified) order and does
bracket selection on the result: BracketingNthOrderBrentSolver.
JIRA: MATH-635
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
(with props)
Modified:
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java?rev=1152276&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
Fri Jul 29 15:44:21 2011
@@ -0,0 +1,397 @@
+/*
+ * 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.math.analysis.solvers;
+
+
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.exception.MathInternalError;
+import org.apache.commons.math.exception.NoBracketingException;
+import org.apache.commons.math.exception.NumberIsTooSmallException;
+import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
+
+/**
+ * This class implements a modification of the <a
+ * href="http://mathworld.wolfram.com/BrentsMethod.html";> Brent algorithm</a>.
+ * <p>
+ * The changes with respect to the original Brent algorithm are:
+ * <ul>
+ * <li>the returned value is chosen in the current interval according
+ * to user specified {@link AllowedSolutions},</li>
+ * <li>the maximal order for the invert polynomial root search is
+ * user-specified instead of being invert quadratic only</li>
+ * </ul>
+ * </p>
+ * The given interval must bracket the root.
+ *
+ * @version $Id$
+ */
+public class BracketingNthOrderBrentSolver
+ extends AbstractUnivariateRealSolver
+ implements BracketedUnivariateRealSolver<UnivariateRealFunction> {
+
+ /** Default absolute accuracy. */
+ private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
+
+ /** Default maximal order. */
+ private static final int DEFAULT_MAXIMAL_ORDER = 5;
+
+ /** Maximal aging triggering an attempt to balance the bracketing
interval. */
+ private static final int MAXIMAL_AGING = 2;
+
+ /** Reduction factor for attempts to balance the bracketing interval. */
+ private static final double REDUCTION_FACTOR = 1.0 / 16.0;
+
+ /** Maximal order. */
+ private final int maximalOrder;
+
+ /** The kinds of solutions that the algorithm may accept. */
+ private AllowedSolutions allowed;
+
+ /**
+ * Construct a solver with default accuracy and maximal order (1e-6 and 5
respectively)
+ */
+ public BracketingNthOrderBrentSolver() {
+ this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
+ }
+
+ /**
+ * Construct a solver.
+ *
+ * @param absoluteAccuracy Absolute accuracy.
+ * @param maximalOrder maximal order.
+ * @exception NumberIsTooSmallException if maximal order is lower than 2
+ */
+ public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
+ final int maximalOrder)
+ throws NumberIsTooSmallException {
+ super(absoluteAccuracy);
+ if (maximalOrder < 2) {
+ throw new NumberIsTooSmallException(maximalOrder, 2, true);
+ }
+ this.maximalOrder = maximalOrder;
+ this.allowed = AllowedSolutions.ANY_SIDE;
+ }
+
+ /**
+ * Construct a solver.
+ *
+ * @param relativeAccuracy Relative accuracy.
+ * @param absoluteAccuracy Absolute accuracy.
+ * @param maximalOrder maximal order.
+ * @exception NumberIsTooSmallException if maximal order is lower than 2
+ */
+ public BracketingNthOrderBrentSolver(final double relativeAccuracy,
+ final double absoluteAccuracy,
+ final int maximalOrder)
+ throws NumberIsTooSmallException {
+ super(relativeAccuracy, absoluteAccuracy);
+ if (maximalOrder < 2) {
+ throw new NumberIsTooSmallException(maximalOrder, 2, true);
+ }
+ this.maximalOrder = maximalOrder;
+ this.allowed = AllowedSolutions.ANY_SIDE;
+ }
+
+ /**
+ * Construct a solver.
+ *
+ * @param relativeAccuracy Relative accuracy.
+ * @param absoluteAccuracy Absolute accuracy.
+ * @param functionValueAccuracy Function value accuracy.
+ * @param maximalOrder maximal order.
+ * @exception NumberIsTooSmallException if maximal order is lower than 2
+ */
+ public BracketingNthOrderBrentSolver(final double relativeAccuracy,
+ final double absoluteAccuracy,
+ final double functionValueAccuracy,
+ final int maximalOrder)
+ throws NumberIsTooSmallException {
+ super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+ if (maximalOrder < 2) {
+ throw new NumberIsTooSmallException(maximalOrder, 2, true);
+ }
+ this.maximalOrder = maximalOrder;
+ this.allowed = AllowedSolutions.ANY_SIDE;
+ }
+
+ /** Get the maximal order.
+ * @return maximal order
+ */
+ public int getMaximalOrder() {
+ return maximalOrder;
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ protected double doSolve() {
+
+ // prepare arrays with the first points
+ final double[] x = new double[maximalOrder + 1];
+ final double[] y = new double[maximalOrder + 1];
+ x[0] = getMin();
+ x[1] = getStartValue();
+ x[2] = getMax();
+ verifySequence(x[0], x[1], x[2]);
+
+ // evaluate initial guess
+ y[1] = computeObjectiveValue(x[1]);
+ if (MathUtils.equals(y[1], 0.0, 1)) {
+ // return the initial guess if it is a perfect root.
+ return x[1];
+ }
+
+ // evaluate first endpoint
+ y[0] = computeObjectiveValue(x[0]);
+ if (MathUtils.equals(y[0], 0.0, 1)) {
+ // return the first endpoint if it is a perfect root.
+ return x[0];
+ }
+
+ int nbPoints;
+ int signChangeIndex;
+ if (y[0] * y[1] < 0) {
+
+ // reduce interval if it brackets the root
+ nbPoints = 2;
+ signChangeIndex = 1;
+
+ } else {
+
+ // evaluate second endpoint
+ y[2] = computeObjectiveValue(x[2]);
+ if (MathUtils.equals(y[2], 0.0, 1)) {
+ // return the second endpoint if it is a perfect root.
+ return x[2];
+ }
+
+ if (y[1] * y[2] < 0) {
+ // use all computed point as a start sampling array for solving
+ nbPoints = 3;
+ signChangeIndex = 2;
+ } else {
+ throw new NoBracketingException(x[0], x[2], y[0], y[2]);
+ }
+
+ }
+
+ // prepare a work array for inverse polynomial interpolation
+ final double[] tmpX = new double[x.length];
+
+ // current tightest bracketing of the root
+ double xA = x[signChangeIndex - 1];
+ double yA = y[signChangeIndex - 1];
+ double absYA = FastMath.abs(yA);
+ int agingA = 0;
+ double xB = x[signChangeIndex];
+ double yB = y[signChangeIndex];
+ double absYB = FastMath.abs(yB);
+ int agingB = 0;
+
+ // search loop
+ while (true) {
+
+ // check convergence of bracketing interval
+ final double xTol = getAbsoluteAccuracy() +
+ getRelativeAccuracy() *
FastMath.max(FastMath.abs(xA), FastMath.abs(xB));
+ if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) <
getFunctionValueAccuracy())) {
+ switch (allowed) {
+ case ANY_SIDE :
+ return absYA < absYB ? xA : xB;
+ case LEFT_SIDE :
+ return xA;
+ case RIGHT_SIDE :
+ return xB;
+ case BELOW_SIDE :
+ return (yA <= 0) ? xA : xB;
+ case ABOVE_SIDE :
+ return (yA < 0) ? xB : xA;
+ default :
+ // this should never happen
+ throw new MathInternalError(null);
+ }
+ }
+
+ // target for the next evaluation point
+ double targetY;
+ if (agingA >= MAXIMAL_AGING) {
+ // we keep updating the high bracket, try to compensate this
+ targetY = -REDUCTION_FACTOR * yB;
+ } else if (agingB >= MAXIMAL_AGING) {
+ // we keep updating the low bracket, try to compensate this
+ targetY = -REDUCTION_FACTOR * yA;
+ } else {
+ // bracketing is balanced, try to find the root itself
+ targetY = 0;
+ }
+
+ // make a few attempts to guess a root,
+ double nextX;
+ int start = 0;
+ int end = nbPoints;
+ do {
+
+ // guess a value for current target, using inverse polynomial
interpolation
+ System.arraycopy(x, start, tmpX, start, end - start);
+ nextX = guessX(targetY, tmpX, y, start, end);
+
+ if (!((nextX > xA) && (nextX < xB))) {
+ // the guessed root is not strictly inside of the tightest
bracketing interval
+
+ // the guessed root is either not strictly inside the
interval or it
+ // is a NaN (which occurs when some sampling points share
the same y)
+ // we try again with a lower interpolation order
+ if (signChangeIndex - start >= end - signChangeIndex) {
+ // we have more points before the sign change, drop
the lowest point
+ ++start;
+ } else {
+ // we have more points after sign change, drop the
highest point
+ --end;
+ }
+
+ // we need to do one more attempt
+ nextX = Double.NaN;
+
+ }
+
+ } while (Double.isNaN(nextX) && (end - start > 1));
+
+ if (Double.isNaN(nextX)) {
+ // fall back to bisection
+ nextX = xA + 0.5 * (xB - xA);
+ start = signChangeIndex - 1;
+ end = signChangeIndex;
+ }
+
+ // evaluate the function at the guessed root
+ final double nextY = computeObjectiveValue(nextX);
+ if (MathUtils.equals(nextY, 0.0, 1)) {
+ // we have found an exact root, since it is not an
approximation
+ // we don't need to bother about the allowed solutions setting
+ return nextX;
+ }
+
+ if ((nbPoints > 2) && (end - start != nbPoints)) {
+
+ // we have been forced to ignore some points to keep
bracketing,
+ // they are probably too far from the root, drop them from now
on
+ nbPoints = end - start;
+ System.arraycopy(x, start, x, 0, nbPoints);
+ System.arraycopy(y, start, y, 0, nbPoints);
+ signChangeIndex -= start;
+
+ } else if (nbPoints == x.length) {
+
+ // we have to drop one point in order to insert the new one
+ nbPoints--;
+
+ // keep the tightest bracketing interval as centered as
possible
+ if (signChangeIndex >= (x.length + 1) / 2) {
+ // we drop the lowest point, we have to shift the arrays
and the index
+ System.arraycopy(x, 1, x, 0, nbPoints);
+ System.arraycopy(y, 1, y, 0, nbPoints);
+ --signChangeIndex;
+ }
+
+ }
+
+ // insert the last computed point
+ //(by construction, we know it lies inside the tightest bracketing
interval)
+ System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1,
nbPoints - signChangeIndex);
+ x[signChangeIndex] = nextX;
+ System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1,
nbPoints - signChangeIndex);
+ y[signChangeIndex] = nextY;
+ ++nbPoints;
+
+ // update the bracketing interval
+ if (nextY * yA <= 0) {
+ // the sign change occurs before the inserted point
+ xB = nextX;
+ yB = nextY;
+ absYB = FastMath.abs(yB);
+ ++agingA;
+ agingB = 0;
+ } else {
+ // the sign change occurs after the inserted point
+ xA = nextX;
+ yA = nextY;
+ absYA = FastMath.abs(yA);
+ agingA = 0;
+ ++agingB;
+
+ // update the sign change index
+ signChangeIndex++;
+
+ }
+
+ }
+
+ }
+
+ /** Guess an x value by n<sup>th</sup> order inverse polynomial
interpolation.
+ * <p>
+ * The x value is guessed by evaluating polynomial Q(y) at y = targetY,
where Q
+ * is built such that for all considered points (x<sub>i</sub>,
y<sub>i</sub>),
+ * Q(y<sub>i</sub>) = x<sub>i</sub>.
+ * </p>
+ * @param targetY target value for y
+ * @param x reference points abscissas for interpolation,
+ * note that this array <em>is</em> modified during computation
+ * @param y reference points ordinates for interpolation
+ * @param start start index of the points to consider (inclusive)
+ * @param end end index of the points to consider (exclusive)
+ * @return guessed root (will be a NaN if two points share the same y)
+ */
+ private double guessX(final double targetY, final double[] x, final
double[] y,
+ final int start, final int end) {
+
+ // compute Q Newton coefficients by divided differences
+ for (int i = start; i < end - 1; ++i) {
+ final int delta = i + 1 - start;
+ for (int j = end - 1; j > i; --j) {
+ x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]);
+ }
+ }
+
+ // evaluate Q(targetY)
+ double x0 = 0;
+ for (int j = end - 1; j >= start; --j) {
+ x0 = x[j] + x0 * (targetY - y[j]);
+ }
+
+ return x0;
+
+ }
+
+ /** {@inheritDoc} */
+ public double solve(int maxEval, UnivariateRealFunction f, double min,
+ double max, AllowedSolutions allowedSolutions) {
+ this.allowed = allowedSolutions;
+ return super.solve(maxEval, f, min, max);
+ }
+
+ /** {@inheritDoc} */
+ public double solve(int maxEval, UnivariateRealFunction f, double min,
+ double max, double startValue,
+ AllowedSolutions allowedSolutions) {
+ this.allowed = allowedSolutions;
+ return super.solve(maxEval, f, min, max, startValue);
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1152276&r1=1152275&r2=1152276&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Jul 29 15:44:21 2011
@@ -52,6 +52,10 @@ The <action> type attribute can be add,u
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="3.0" date="TBD" description="TBD">
+ <action dev="luc" type="add" issue="MATH-635" >
+ Added a Brent-like solver that has higher (user specified) order and
does
+ bracket selection on the result: BracketingNthOrderBrentSolver.
+ </action>
<action dev="luc" type="add" >
Added a few shortcut methods and predicates to Dfp (abs, isZero,
negativeOrNull, strictlyNegative, positiveOrNull, strictlyPositive).
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml?rev=1152276&r1=1152275&r2=1152276&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml Fri Jul 29
15:44:21 2011
@@ -74,6 +74,12 @@
<td>no</td>
</tr>
<tr>
+ <td><a
href="../apidocs/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.html">bracketing
n<sup>th</sup> order Brent</a></td>
+ <td>variable order, guaranteed</td>
+ <td>yes</td>
+ <td>yes</td>
+ </tr>
+ <tr>
<td><a
href="../apidocs/org/apache/commons/math/analysis/solvers/IllinoisSolver.html">Illinois
Method</a></td>
<td><a
href="../apidocs/org/apache/commons/math/analysis/UnivariateRealFunction.html">univariate
real-valued functions</a></td>
<td>super-linear, guaranteed</td>
@@ -206,7 +212,8 @@
<source>UnivariateRealFunction function = // some user defined
function object
final double relativeAccuracy = 1.0e-12;
final double absoluteAccuracy = 1.0e-8;
-UnivariateRealSolver solver = new PegasusSolver(relativeAccuracy,
absoluteAccuracy);
+final int maxOrder = 5;
+UnivariateRealSolver solver = new
BracketingNthOrderBrentSolver(relativeAccuracy, absoluteAccuracy, maxOrder);
double c = solver.solve(100, function, 1.0, 5.0,
AllowedSolutions.LEFT_SIDE);</source>
<p>
Force bracketing, by refining a base solution found by a
non-bracketing solver:
@@ -221,9 +228,9 @@ double c = UnivariateRealSolverUtils.for
baseRoot, 1.0, 5.0,
AllowedSolutions.LEFT_SIDE);
</source>
<p>
- The <code>BrentSolve</code> uses the Brent-Dekker algorithm which is
- fast and robust. This algorithm is recommended for most users. If
there
- are multiple roots in the interval, or there is a large domain of
indeterminacy, the
+ The <code>BrentSolver</code> uses the Brent-Dekker algorithm which is
+ fast and robust. If there are multiple roots in the interval,
+ or there is a large domain of indeterminacy, the
algorithm will converge to a random root in the interval without
indication that there are problems. Interestingly, the examined text
book implementations all disagree in details of the convergence
@@ -233,6 +240,15 @@ double c = UnivariateRealSolverUtils.for
algorithm.
</p>
<p>
+ The <code>BracketingNthOrderBrentSolver</code> uses an extension of
the
+ Brent-Dekker algorithm which uses inverse n<sup>th</sup> order
polynomial
+ interpolation instead of inverse quadratic interpolation, and which
allows
+ selection of the side of the convergence interval for result
bracketing.
+ This is now the recommended algorithm for most users since it has the
+ largest order, doesn't require derivatives, has guaranteed
convergence
+ and allows result bracket selection.
+ </p>
+ <p>
The <code>SecantSolver</code> uses a straightforward secant
algorithm which does not bracket the search and therefore does not
guarantee convergence. It may be faster than Brent on some
well-behaved
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java?rev=1152276&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
Fri Jul 29 15:44:21 2011
@@ -0,0 +1,181 @@
+/*
+ * 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.math.analysis.solvers;
+
+import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
+import org.apache.commons.math.analysis.QuinticFunction;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.exception.NumberIsTooSmallException;
+import org.apache.commons.math.exception.TooManyEvaluationsException;
+import org.apache.commons.math.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test case for {@link BracketingNthOrderBrentSolver bracketing
n<sup>th</sup> order Brent} solver.
+ *
+ * @version $Id$
+ */
+public final class BracketingNthOrderBrentSolverTest extends
BaseSecantSolverAbstractTest {
+ /** {@inheritDoc} */
+ protected UnivariateRealSolver getSolver() {
+ return new BracketingNthOrderBrentSolver();
+ }
+
+ /** {@inheritDoc} */
+ protected int[] getQuinticEvalCounts() {
+ return new int[] {1, 3, 8, 1, 9, 4, 8, 1, 12, 1, 14};
+ }
+
+ @Test(expected=NumberIsTooSmallException.class)
+ public void testInsufficientOrder1() {
+ new BracketingNthOrderBrentSolver(1.0e-10, 1);
+ }
+
+ @Test(expected=NumberIsTooSmallException.class)
+ public void testInsufficientOrder2() {
+ new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1);
+ }
+
+ @Test(expected=NumberIsTooSmallException.class)
+ public void testInsufficientOrder3() {
+ new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 1);
+ }
+
+ @Test
+ public void testConstructorsOK() {
+ Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10,
2).getMaximalOrder());
+ Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10,
1.0e-10, 2).getMaximalOrder());
+ Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10,
1.0e-10, 1.0e-10, 2).getMaximalOrder());
+ }
+
+ @Test
+ public void testConvergenceOnFunctionAccuracy() {
+ BracketingNthOrderBrentSolver solver =
+ new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-10, 0.001, 3);
+ QuinticFunction f = new QuinticFunction();
+ double result = solver.solve(20, f, 0.2, 0.9, 0.4,
AllowedSolutions.BELOW_SIDE);
+ Assert.assertEquals(0, f.value(result),
solver.getFunctionValueAccuracy());
+ Assert.assertTrue(f.value(result) <= 0);
+ Assert.assertTrue(result - 0.5 > solver.getAbsoluteAccuracy());
+ result = solver.solve(20, f, -0.9, -0.2, -0.4,
AllowedSolutions.ABOVE_SIDE);
+ Assert.assertEquals(0, f.value(result),
solver.getFunctionValueAccuracy());
+ Assert.assertTrue(f.value(result) >= 0);
+ Assert.assertTrue(result + 0.5 < -solver.getAbsoluteAccuracy());
+ }
+
+ @Test
+ public void testFasterThanNewton() {
+ // the following test functions come from Beny Neta's paper:
+ // "Several New Methods for solving Equations"
+ // intern J. Computer Math Vol 23 pp 265-282
+ // available here:
http://www.math.nps.navy.mil/~bneta/SeveralNewMethods.PDF
+ // the reference roots have been computed by the Dfp solver to more
than
+ // 80 digits and checked with emacs (only the first 20 digits are
reproduced here)
+ compare(new TestFunction(0.0, -2, 2) {
+ public double value(double x) { return FastMath.sin(x) - 0.5
* x; }
+ public double derivative(double x) { return FastMath.cos(x) - 0.5;
}
+ });
+ compare(new TestFunction(6.3087771299726890947, -5, 10) {
+ public double value(double x) { return FastMath.pow(x, 5) + x
- 10000; }
+ public double derivative(double x) { return 5 * FastMath.pow(x, 4)
+ 1; }
+ });
+ compare(new TestFunction(9.6335955628326951924, 0.001, 10) {
+ public double value(double x) { return FastMath.sqrt(x) - 1 /
x - 3; }
+ public double derivative(double x) { return 0.5 / FastMath.sqrt(x)
+ 1 / (x * x); }
+ });
+ compare(new TestFunction(2.8424389537844470678, -5, 5) {
+ public double value(double x) { return FastMath.exp(x) + x -
20; }
+ public double derivative(double x) { return FastMath.exp(x) + 1; }
+ });
+ compare(new TestFunction(8.3094326942315717953, 0.001, 10) {
+ public double value(double x) { return FastMath.log(x) +
FastMath.sqrt(x) - 5; }
+ public double derivative(double x) { return 1 / x + 0.5 /
FastMath.sqrt(x); }
+ });
+ compare(new TestFunction(1.4655712318767680266, -0.5, 1.5) {
+ public double value(double x) { return (x - 1) * x * x - 1; }
+ public double derivative(double x) { return (3 * x - 2) * x; }
+ });
+
+ }
+
+ private void compare(TestFunction f) {
+ compare(f, f.getRoot(), f.getMin(), f.getMax());
+ }
+
+ private void compare(DifferentiableUnivariateRealFunction f,
+ double root, double min, double max) {
+ NewtonSolver newton = new NewtonSolver(1.0e-12);
+ BracketingNthOrderBrentSolver bracketing =
+ new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-12, 1.0e-18,
5);
+ double resultN;
+ try {
+ resultN = newton.solve(100, f, min, max);
+ } catch (TooManyEvaluationsException tmee) {
+ resultN = Double.NaN;
+ }
+ double resultB;
+ try {
+ resultB = bracketing.solve(100, f, min, max);
+ } catch (TooManyEvaluationsException tmee) {
+ resultB = Double.NaN;
+ }
+ Assert.assertEquals(root, resultN, newton.getAbsoluteAccuracy());
+ Assert.assertEquals(root, resultB, bracketing.getAbsoluteAccuracy());
+ Assert.assertTrue(bracketing.getEvaluations() <
newton.getEvaluations());
+ }
+
+ private static abstract class TestFunction implements
DifferentiableUnivariateRealFunction {
+
+ private final double root;
+ private final double min;
+ private final double max;
+
+ protected TestFunction(final double root, final double min, final
double max) {
+ this.root = root;
+ this.min = min;
+ this.max = max;
+ }
+
+ public double getRoot() {
+ return root;
+ }
+
+ public double getMin() {
+ return min;
+ }
+
+ public double getMax() {
+ return max;
+ }
+
+ public abstract double value(double x);
+
+ public abstract double derivative(double x);
+
+ public UnivariateRealFunction derivative() {
+ return new UnivariateRealFunction() {
+ public double value(double x) {
+ return derivative(x);
+ }
+ };
+ }
+
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
