Added a Brent-like solver that has higher (user specified) order and does
bracket selection on the result: BracketingNthOrderBrentSolver. JIRA: MATH-635 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1152276 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
21d9fa79fb
commit
c68662e257
|
|
@ -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);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -52,6 +52,10 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
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).
|
||||
|
|
|
|||
|
|
@ -73,6 +73,12 @@
|
|||
<td>yes</td>
|
||||
<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>
|
||||
|
|
@ -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.forceSide(100, function,
|
|||
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
|
||||
|
|
@ -232,6 +239,15 @@ double c = UnivariateRealSolverUtils.forceSide(100, function,
|
|||
get exactly the same root values as for other implementations of this
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue