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

Added and used a specialized convergence exception for exceeded iteration counts 

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@506585 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2007-02-12 19:17:41 +00:00
parent 32ea2a389a
commit ee71801e77
14 changed files with 903 additions and 817 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
src/java/org/apache/commons/math/MaxIterationsExceededException.java Normal file
View File

@ -0,0 +1,44 @@
package org.apache.commons.math;
import org.apache.commons.math.ConvergenceException;
public class MaxIterationsExceededException extends ConvergenceException {
/** Serializable version identifier. */
private static final long serialVersionUID = -2154780004193976271L;
/** Maximal number of iterations allowed. */
private int maxIterations;
/**
* Constructs an exception with specified formatted detail message.
* Message formatting is delegated to {@link java.text.MessageFormat}.
* @param maxIterations maximal number of iterations allowed
*/
public MaxIterationsExceededException(int maxIterations) {
super("Maximal number of iterations ({0}) exceeded",
new Object[] { new Integer(maxIterations) });
this.maxIterations = maxIterations;
}
/**
* Constructs an exception with specified formatted detail message.
* Message formatting is delegated to {@link java.text.MessageFormat}.
* @param argument the failing function argument
* @param pattern format specifier
* @param arguments format arguments
*/
public MaxIterationsExceededException(int maxIterations,
String pattern, Object[] arguments) {
super(pattern, arguments);
this.maxIterations = maxIterations;
}
/** Get the maximal number of iterations allowed.
* @return maximal number of iterations allowed
*/
public int getMaxIterations() {
return maxIterations;
}
}

17
src/java/org/apache/commons/math/analysis/BisectionSolver.java
View File

@ -17,7 +17,7 @@
package org.apache.commons.math.analysis;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements the <a href="http://mathworld.wolfram.com/Bisection.html">
@ -30,8 +30,8 @@ import org.apache.commons.math.ConvergenceException;
public class BisectionSolver extends UnivariateRealSolverImpl {
/** Serializable version identifier */
private static final long serialVersionUID = 7137520585963699578L;
private static final long serialVersionUID = 4963578633786538912L;
/**
* Construct a solver for the given function.
*
@ -48,13 +48,13 @@ public class BisectionSolver extends UnivariateRealSolverImpl {
* @param max the upper bound for the interval.
* @param initial the start value to use (ignored).
* @return the value where the function is zero
* @throws ConvergenceException the maximum iteration count is exceeded
* @throws MaxIterationsExceededException the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating
* the function
* @throws IllegalArgumentException if min is not less than max
*/
public double solve(double min, double max, double initial)
throws ConvergenceException, FunctionEvaluationException {
throws MaxIterationsExceededException, FunctionEvaluationException {
return solve(min, max);
}
@ -65,12 +65,12 @@ public class BisectionSolver extends UnivariateRealSolverImpl {
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded.
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded.
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min is not less than max
*/
public double solve(double min, double max) throws ConvergenceException,
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
clearResult();
@ -101,7 +101,6 @@ public class BisectionSolver extends UnivariateRealSolverImpl {
++i;
}
throw new ConvergenceException
("Maximum number of iterations exceeded: " + maximalIterationCount);
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

14
src/java/org/apache/commons/math/analysis/BrentSolver.java
View File

@ -17,8 +17,8 @@
package org.apache.commons.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
@ -31,7 +31,7 @@ import org.apache.commons.math.FunctionEvaluationException;
public class BrentSolver extends UnivariateRealSolverImpl {
/** Serializable version identifier */
private static final long serialVersionUID = 3350616277306882875L;
private static final long serialVersionUID = -2136672307739067002L;
/**
* Construct a solver for the given function.
@ -52,13 +52,13 @@ public class BrentSolver extends UnivariateRealSolverImpl {
* @param max the upper bound for the interval.
* @param initial the start value to use (ignored).
* @return the value where the function is zero
* @throws ConvergenceException the maximum iteration count is exceeded
* @throws MaxIterationsExceededException the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating
* the function
* @throws IllegalArgumentException if initial is not between min and max
*/
public double solve(double min, double max, double initial)
throws ConvergenceException, FunctionEvaluationException {
throws MaxIterationsExceededException, FunctionEvaluationException {
return solve(min, max);
}
@ -73,13 +73,13 @@ public class BrentSolver extends UnivariateRealSolverImpl {
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
*/
public double solve(double min, double max) throws ConvergenceException,
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
clearResult();
@ -189,6 +189,6 @@ public class BrentSolver extends UnivariateRealSolverImpl {
}
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

653
src/java/org/apache/commons/math/analysis/LaguerreSolver.java Executable file → Normal file
View File

@ -1,326 +1,327 @@
/*
* 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;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.complex.Complex;
import org.apache.commons.math.complex.ComplexUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
* Laguerre's Method</a> for root finding of real coefficient polynomials.
* For reference, see <b>A First Course in Numerical Analysis</b>,
* ISBN 048641454X, chapter 8.
* <p>
* Laguerre's method is global in the sense that it can start with any initial
* approximation and be able to solve all roots from that point.
*
* @version $Revision$ $Date$
*/
public class LaguerreSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
private static final long serialVersionUID = 5287689975005870178L;
/** polynomial function to solve */
private PolynomialFunction p;
/**
* Construct a solver for the given function.
*
* @param f function to solve
* @throws IllegalArgumentException if function is not polynomial
*/
public LaguerreSolver(UnivariateRealFunction f) throws
IllegalArgumentException {
super(f, 100, 1E-6);
if (f instanceof PolynomialFunction) {
p = (PolynomialFunction)f;
} else {
throw new IllegalArgumentException("Function is not polynomial.");
}
}
/**
* Returns a copy of the polynomial function.
*
* @return a fresh copy of the polynomial function
*/
public PolynomialFunction getPolynomialFunction() {
return new PolynomialFunction(p.getCoefficients());
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
ConvergenceException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (p.value(min) == 0.0) { return min; }
if (p.value(max) == 0.0) { return max; }
if (p.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, p);
verifySequence(min, initial, max);
if (isBracketing(min, initial, p)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a real root in the given interval.
* <p>
* Despite the bracketing condition, the root returned by solve(Complex[],
* Complex) may not be a real zero inside [min, max]. For example,
* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
* another initial value, or, as we did here, call solveAll() to obtain
* all roots and pick up the one that we're looking for.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws ConvergenceException,
FunctionEvaluationException {
// check for zeros before verifying bracketing
if (p.value(min) == 0.0) { return min; }
if (p.value(max) == 0.0) { return max; }
verifyBracketing(min, max, p);
double coefficients[] = p.getCoefficients();
Complex c[] = new Complex[coefficients.length];
for (int i = 0; i < coefficients.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
}
Complex initial = new Complex(0.5 * (min + max), 0.0);
Complex z = solve(c, initial);
if (isRootOK(min, max, z)) {
setResult(z.getReal(), iterationCount);
return result;
}
// solve all roots and select the one we're seeking
Complex[] root = solveAll(c, initial);
for (int i = 0; i < root.length; i++) {
if (isRootOK(min, max, root[i])) {
setResult(root[i].getReal(), iterationCount);
return result;
}
}
// should never happen
throw new ConvergenceException("Convergence failed.");
}
/**
* Returns true iff the given complex root is actually a real zero
* in the given interval, within the solver tolerance level.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param z the complex root
* @return true iff z is the sought-after real zero
*/
protected boolean isRootOK(double min, double max, Complex z) {
double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
return (isSequence(min, z.getReal(), max)) &&
(Math.abs(z.getImaginary()) <= tolerance ||
z.abs() <= functionValueAccuracy);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(double coefficients[], double initial) throws
ConvergenceException, FunctionEvaluationException {
Complex c[] = new Complex[coefficients.length];
Complex z = new Complex(initial, 0.0);
for (int i = 0; i < c.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
}
return solveAll(c, z);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(Complex coefficients[], Complex initial) throws
ConvergenceException, FunctionEvaluationException {
int n = coefficients.length - 1;
int iterationCount = 0;
if (n < 1) {
throw new IllegalArgumentException
("Polynomial degree must be positive: degree=" + n);
}
Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial
for (int i = 0; i <= n; i++) {
c[i] = coefficients[i];
}
// solve individual root successively
Complex root[] = new Complex[n];
for (int i = 0; i < n; i++) {
Complex subarray[] = new Complex[n-i+1];
System.arraycopy(c, 0, subarray, 0, subarray.length);
root[i] = solve(subarray, initial);
// polynomial deflation using synthetic division
Complex newc = c[n-i];
Complex oldc = null;
for (int j = n-i-1; j >= 0; j--) {
oldc = c[j];
c[j] = newc;
newc = oldc.add(newc.multiply(root[i]));
}
iterationCount += this.iterationCount;
}
resultComputed = true;
this.iterationCount = iterationCount;
return root;
}
/**
* Find a complex root for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex solve(Complex coefficients[], Complex initial) throws
ConvergenceException, FunctionEvaluationException {
int n = coefficients.length - 1;
if (n < 1) {
throw new IllegalArgumentException
("Polynomial degree must be positive: degree=" + n);
}
Complex N = new Complex((double)n, 0.0);
Complex N1 = new Complex((double)(n-1), 0.0);
int i = 1;
Complex pv = null;
Complex dv = null;
Complex d2v = null;
Complex G = null;
Complex G2 = null;
Complex H = null;
Complex delta = null;
Complex denominator = null;
Complex z = initial;
Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
while (i <= maximalIterationCount) {
// Compute pv (polynomial value), dv (derivative value), and
// d2v (second derivative value) simultaneously.
pv = coefficients[n];
dv = Complex.ZERO;
d2v = Complex.ZERO;
for (int j = n-1; j >= 0; j--) {
d2v = dv.add(z.multiply(d2v));
dv = pv.add(z.multiply(dv));
pv = coefficients[j].add(z.multiply(pv));
}
d2v = d2v.multiply(new Complex(2.0, 0.0));
// check for convergence
double tolerance = Math.max(relativeAccuracy * z.abs(),
absoluteAccuracy);
if ((z.subtract(oldz)).abs() <= tolerance) {
resultComputed = true;
iterationCount = i;
return z;
}
if (pv.abs() <= functionValueAccuracy) {
resultComputed = true;
iterationCount = i;
return z;
}
// now pv != 0, calculate the new approximation
G = dv.divide(pv);
G2 = G.multiply(G);
H = G2.subtract(d2v.divide(pv));
delta = N1.multiply((N.multiply(H)).subtract(G2));
// choose a denominator larger in magnitude
Complex dplus = G.add(ComplexUtils.sqrt(delta));
Complex dminus = G.subtract(ComplexUtils.sqrt(delta));
denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
// Perturb z if denominator is zero, for instance,
// p(x) = x^3 + 1, z = 0.
if (denominator.equals(new Complex(0.0, 0.0))) {
z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
oldz = new Complex(Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY);
} else {
oldz = z;
z = z.subtract(N.divide(denominator));
}
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
}
/*
* 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;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.complex.Complex;
import org.apache.commons.math.complex.ComplexUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
* Laguerre's Method</a> for root finding of real coefficient polynomials.
* For reference, see <b>A First Course in Numerical Analysis</b>,
* ISBN 048641454X, chapter 8.
* <p>
* Laguerre's method is global in the sense that it can start with any initial
* approximation and be able to solve all roots from that point.
*
* @version $Revision$ $Date$
*/
public class LaguerreSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
private static final long serialVersionUID = -3775334783473775723L;
/** polynomial function to solve */
private PolynomialFunction p;
/**
* Construct a solver for the given function.
*
* @param f function to solve
* @throws IllegalArgumentException if function is not polynomial
*/
public LaguerreSolver(UnivariateRealFunction f) throws
IllegalArgumentException {
super(f, 100, 1E-6);
if (f instanceof PolynomialFunction) {
p = (PolynomialFunction)f;
} else {
throw new IllegalArgumentException("Function is not polynomial.");
}
}
/**
* Returns a copy of the polynomial function.
*
* @return a fresh copy of the polynomial function
*/
public PolynomialFunction getPolynomialFunction() {
return new PolynomialFunction(p.getCoefficients());
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
ConvergenceException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (p.value(min) == 0.0) { return min; }
if (p.value(max) == 0.0) { return max; }
if (p.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, p);
verifySequence(min, initial, max);
if (isBracketing(min, initial, p)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a real root in the given interval.
* <p>
* Despite the bracketing condition, the root returned by solve(Complex[],
* Complex) may not be a real zero inside [min, max]. For example,
* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
* another initial value, or, as we did here, call solveAll() to obtain
* all roots and pick up the one that we're looking for.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws ConvergenceException,
FunctionEvaluationException {
// check for zeros before verifying bracketing
if (p.value(min) == 0.0) { return min; }
if (p.value(max) == 0.0) { return max; }
verifyBracketing(min, max, p);
double coefficients[] = p.getCoefficients();
Complex c[] = new Complex[coefficients.length];
for (int i = 0; i < coefficients.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
}
Complex initial = new Complex(0.5 * (min + max), 0.0);
Complex z = solve(c, initial);
if (isRootOK(min, max, z)) {
setResult(z.getReal(), iterationCount);
return result;
}
// solve all roots and select the one we're seeking
Complex[] root = solveAll(c, initial);
for (int i = 0; i < root.length; i++) {
if (isRootOK(min, max, root[i])) {
setResult(root[i].getReal(), iterationCount);
return result;
}
}
// should never happen
throw new ConvergenceException();
}
/**
* Returns true iff the given complex root is actually a real zero
* in the given interval, within the solver tolerance level.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param z the complex root
* @return true iff z is the sought-after real zero
*/
protected boolean isRootOK(double min, double max, Complex z) {
double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
return (isSequence(min, z.getReal(), max)) &&
(Math.abs(z.getImaginary()) <= tolerance ||
z.abs() <= functionValueAccuracy);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(double coefficients[], double initial) throws
ConvergenceException, FunctionEvaluationException {
Complex c[] = new Complex[coefficients.length];
Complex z = new Complex(initial, 0.0);
for (int i = 0; i < c.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
}
return solveAll(c, z);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(Complex coefficients[], Complex initial) throws
MaxIterationsExceededException, FunctionEvaluationException {
int n = coefficients.length - 1;
int iterationCount = 0;
if (n < 1) {
throw new IllegalArgumentException
("Polynomial degree must be positive: degree=" + n);
}
Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial
for (int i = 0; i <= n; i++) {
c[i] = coefficients[i];
}
// solve individual root successively
Complex root[] = new Complex[n];
for (int i = 0; i < n; i++) {
Complex subarray[] = new Complex[n-i+1];
System.arraycopy(c, 0, subarray, 0, subarray.length);
root[i] = solve(subarray, initial);
// polynomial deflation using synthetic division
Complex newc = c[n-i];
Complex oldc = null;
for (int j = n-i-1; j >= 0; j--) {
oldc = c[j];
c[j] = newc;
newc = oldc.add(newc.multiply(root[i]));
}
iterationCount += this.iterationCount;
}
resultComputed = true;
this.iterationCount = iterationCount;
return root;
}
/**
* Find a complex root for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex solve(Complex coefficients[], Complex initial) throws
MaxIterationsExceededException, FunctionEvaluationException {
int n = coefficients.length - 1;
if (n < 1) {
throw new IllegalArgumentException
("Polynomial degree must be positive: degree=" + n);
}
Complex N = new Complex((double)n, 0.0);
Complex N1 = new Complex((double)(n-1), 0.0);
int i = 1;
Complex pv = null;
Complex dv = null;
Complex d2v = null;
Complex G = null;
Complex G2 = null;
Complex H = null;
Complex delta = null;
Complex denominator = null;
Complex z = initial;
Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
while (i <= maximalIterationCount) {
// Compute pv (polynomial value), dv (derivative value), and
// d2v (second derivative value) simultaneously.
pv = coefficients[n];
dv = Complex.ZERO;
d2v = Complex.ZERO;
for (int j = n-1; j >= 0; j--) {
d2v = dv.add(z.multiply(d2v));
dv = pv.add(z.multiply(dv));
pv = coefficients[j].add(z.multiply(pv));
}
d2v = d2v.multiply(new Complex(2.0, 0.0));
// check for convergence
double tolerance = Math.max(relativeAccuracy * z.abs(),
absoluteAccuracy);
if ((z.subtract(oldz)).abs() <= tolerance) {
resultComputed = true;
iterationCount = i;
return z;
}
if (pv.abs() <= functionValueAccuracy) {
resultComputed = true;
iterationCount = i;
return z;
}
// now pv != 0, calculate the new approximation
G = dv.divide(pv);
G2 = G.multiply(G);
H = G2.subtract(d2v.divide(pv));
delta = N1.multiply((N.multiply(H)).subtract(G2));
// choose a denominator larger in magnitude
Complex dplus = G.add(ComplexUtils.sqrt(delta));
Complex dminus = G.subtract(ComplexUtils.sqrt(delta));
denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
// Perturb z if denominator is zero, for instance,
// p(x) = x^3 + 1, z = 0.
if (denominator.equals(new Complex(0.0, 0.0))) {
z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
oldz = new Complex(Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY);
} else {
oldz = z;
z = z.subtract(N.divide(denominator));
}
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

556
src/java/org/apache/commons/math/analysis/MullerSolver.java Executable file → Normal file
View File

@ -1,278 +1,278 @@
/*
* 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;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
* Muller's Method</a> for root finding of real univariate functions. For
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
* chapter 3.
* <p>
* Muller's method applies to both real and complex functions, but here we
* restrict ourselves to real functions. Methods solve() and solve2() find
* real zeros, using different ways to bypass complex arithmetics.
*
* @version $Revision$ $Date$
*/
public class MullerSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
static final long serialVersionUID = 2619993603551148137L;
/**
* Construct a solver for the given function.
*
* @param f function to solve
*/
public MullerSolver(UnivariateRealFunction f) {
super(f, 100, 1E-6);
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
ConvergenceException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a real root in the given interval.
* <p>
* Original Muller's method would have function evaluation at complex point.
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
* condition is one way to go: by requiring bracketing in every iteration,
* the newly computed approximation is guaranteed to be real.
* <p>
* Normally Muller's method converges quadratically in the vicinity of a
* zero, however it may be very slow in regions far away from zeros. For
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
* bisection as a safety backup if it performs very poorly.
* <p>
* The formulas here use divided differences directly.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws ConvergenceException,
FunctionEvaluationException {
// [x0, x2] is the bracketing interval in each iteration
// x1 is the last approximation and an interpolation point in (x0, x2)
// x is the new root approximation and new x1 for next round
// d01, d12, d012 are divided differences
double x0, x1, x2, x, oldx, y0, y1, y2, y;
double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
x0 = min; y0 = f.value(x0);
x2 = max; y2 = f.value(x2);
x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
// check for zeros before verifying bracketing
if (y0 == 0.0) { return min; }
if (y2 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// Muller's method employs quadratic interpolation through
// x0, x1, x2 and x is the zero of the interpolating parabola.
// Due to bracketing condition, this parabola must have two
// real roots and we choose one in [x0, x2] to be x.
d01 = (y1 - y0) / (x1 - x0);
d12 = (y2 - y1) / (x2 - x1);
d012 = (d12 - d01) / (x2 - x0);
c1 = d01 + (x1 - x0) * d012;
delta = c1 * c1 - 4 * y1 * d012;
xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
// xplus and xminus are two roots of parabola and at least
// one of them should lie in (x0, x2)
x = isSequence(x0, xplus, x2) ? xplus : xminus;
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// Bisect if convergence is too slow. Bisection would waste
// our calculation of x, hopefully it won't happen often.
boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
(x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
(x == x1);
// prepare the new bracketing interval for next iteration
if (!bisect) {
x0 = x < x1 ? x0 : x1;
y0 = x < x1 ? y0 : y1;
x2 = x > x1 ? x2 : x1;
y2 = x > x1 ? y2 : y1;
x1 = x; y1 = y;
oldx = x;
} else {
double xm = 0.5 * (x0 + x2);
double ym = f.value(xm);
if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
x2 = xm; y2 = ym;
} else {
x0 = xm; y0 = ym;
}
x1 = 0.5 * (x0 + x2);
y1 = f.value(x1);
oldx = Double.POSITIVE_INFINITY;
}
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
/**
* Find a real root in the given interval.
* <p>
* solve2() differs from solve() in the way it avoids complex operations.
* Except for the initial [min, max], solve2() does not require bracketing
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
* number arises in the computation, we simply use its modulus as real
* approximation.
* <p>
* Because the interval may not be bracketing, bisection alternative is
* not applicable here. However in practice our treatment usually works
* well, especially near real zeros where the imaginary part of complex
* approximation is often negligible.
* <p>
* The formulas here do not use divided differences directly.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve2(double min, double max) throws ConvergenceException,
FunctionEvaluationException {
// x2 is the last root approximation
// x is the new approximation and new x2 for next round
// x0 < x1 < x2 does not hold here
double x0, x1, x2, x, oldx, y0, y1, y2, y;
double q, A, B, C, delta, denominator, tolerance;
x0 = min; y0 = f.value(x0);
x1 = max; y1 = f.value(x1);
x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y0 == 0.0) { return min; }
if (y1 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// quadratic interpolation through x0, x1, x2
q = (x2 - x1) / (x1 - x0);
A = q * (y2 - (1 + q) * y1 + q * y0);
B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
C = (1 + q) * y2;
delta = B * B - 4 * A * C;
if (delta >= 0.0) {
// choose a denominator larger in magnitude
double dplus = B + Math.sqrt(delta);
double dminus = B - Math.sqrt(delta);
denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
} else {
// take the modulus of (B +/- Math.sqrt(delta))
denominator = Math.sqrt(B * B - delta);
}
if (denominator != 0) {
x = x2 - 2.0 * C * (x2 - x1) / denominator;
// perturb x if it coincides with x1 or x2
while (x == x1 || x == x2) {
x += absoluteAccuracy;
}
} else {
// extremely rare case, get a random number to skip it
x = min + Math.random() * (max - min);
oldx = Double.POSITIVE_INFINITY;
}
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the next iteration
x0 = x1; y0 = y1;
x1 = x2; y1 = y2;
x2 = x; y2 = y;
oldx = x;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
}
/*
* 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;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
* Muller's Method</a> for root finding of real univariate functions. For
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
* chapter 3.
* <p>
* Muller's method applies to both real and complex functions, but here we
* restrict ourselves to real functions. Methods solve() and solve2() find
* real zeros, using different ways to bypass complex arithmetics.
*
* @version $Revision$ $Date$
*/
public class MullerSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
private static final long serialVersionUID = 6552227503458976920L;
/**
* Construct a solver for the given function.
*
* @param f function to solve
*/
public MullerSolver(UnivariateRealFunction f) {
super(f, 100, 1E-6);
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
MaxIterationsExceededException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a real root in the given interval.
* <p>
* Original Muller's method would have function evaluation at complex point.
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
* condition is one way to go: by requiring bracketing in every iteration,
* the newly computed approximation is guaranteed to be real.
* <p>
* Normally Muller's method converges quadratically in the vicinity of a
* zero, however it may be very slow in regions far away from zeros. For
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
* bisection as a safety backup if it performs very poorly.
* <p>
* The formulas here use divided differences directly.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
// [x0, x2] is the bracketing interval in each iteration
// x1 is the last approximation and an interpolation point in (x0, x2)
// x is the new root approximation and new x1 for next round
// d01, d12, d012 are divided differences
double x0, x1, x2, x, oldx, y0, y1, y2, y;
double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
x0 = min; y0 = f.value(x0);
x2 = max; y2 = f.value(x2);
x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
// check for zeros before verifying bracketing
if (y0 == 0.0) { return min; }
if (y2 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// Muller's method employs quadratic interpolation through
// x0, x1, x2 and x is the zero of the interpolating parabola.
// Due to bracketing condition, this parabola must have two
// real roots and we choose one in [x0, x2] to be x.
d01 = (y1 - y0) / (x1 - x0);
d12 = (y2 - y1) / (x2 - x1);
d012 = (d12 - d01) / (x2 - x0);
c1 = d01 + (x1 - x0) * d012;
delta = c1 * c1 - 4 * y1 * d012;
xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
// xplus and xminus are two roots of parabola and at least
// one of them should lie in (x0, x2)
x = isSequence(x0, xplus, x2) ? xplus : xminus;
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// Bisect if convergence is too slow. Bisection would waste
// our calculation of x, hopefully it won't happen often.
boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
(x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
(x == x1);
// prepare the new bracketing interval for next iteration
if (!bisect) {
x0 = x < x1 ? x0 : x1;
y0 = x < x1 ? y0 : y1;
x2 = x > x1 ? x2 : x1;
y2 = x > x1 ? y2 : y1;
x1 = x; y1 = y;
oldx = x;
} else {
double xm = 0.5 * (x0 + x2);
double ym = f.value(xm);
if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
x2 = xm; y2 = ym;
} else {
x0 = xm; y0 = ym;
}
x1 = 0.5 * (x0 + x2);
y1 = f.value(x1);
oldx = Double.POSITIVE_INFINITY;
}
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**
* Find a real root in the given interval.
* <p>
* solve2() differs from solve() in the way it avoids complex operations.
* Except for the initial [min, max], solve2() does not require bracketing
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
* number arises in the computation, we simply use its modulus as real
* approximation.
* <p>
* Because the interval may not be bracketing, bisection alternative is
* not applicable here. However in practice our treatment usually works
* well, especially near real zeros where the imaginary part of complex
* approximation is often negligible.
* <p>
* The formulas here do not use divided differences directly.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve2(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
// x2 is the last root approximation
// x is the new approximation and new x2 for next round
// x0 < x1 < x2 does not hold here
double x0, x1, x2, x, oldx, y0, y1, y2, y;
double q, A, B, C, delta, denominator, tolerance;
x0 = min; y0 = f.value(x0);
x1 = max; y1 = f.value(x1);
x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y0 == 0.0) { return min; }
if (y1 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// quadratic interpolation through x0, x1, x2
q = (x2 - x1) / (x1 - x0);
A = q * (y2 - (1 + q) * y1 + q * y0);
B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
C = (1 + q) * y2;
delta = B * B - 4 * A * C;
if (delta >= 0.0) {
// choose a denominator larger in magnitude
double dplus = B + Math.sqrt(delta);
double dminus = B - Math.sqrt(delta);
denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
} else {
// take the modulus of (B +/- Math.sqrt(delta))
denominator = Math.sqrt(B * B - delta);
}
if (denominator != 0) {
x = x2 - 2.0 * C * (x2 - x1) / denominator;
// perturb x if it coincides with x1 or x2
while (x == x1 || x == x2) {
x += absoluteAccuracy;
}
} else {
// extremely rare case, get a random number to skip it
x = min + Math.random() * (max - min);
oldx = Double.POSITIVE_INFINITY;
}
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the next iteration
x0 = x1; y0 = y1;
x1 = x2; y1 = y2;
x2 = x; y2 = y;
oldx = x;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

17
src/java/org/apache/commons/math/analysis/NewtonSolver.java
View File

@ -18,8 +18,8 @@
package org.apache.commons.math.analysis;
import java.io.IOException;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html">
@ -32,8 +32,8 @@ import org.apache.commons.math.FunctionEvaluationException;
public class NewtonSolver extends UnivariateRealSolverImpl {
/** Serializable version identifier */
private static final long serialVersionUID = 2606474895443431607L;
private static final long serialVersionUID = 2067325783137941016L;
/** The first derivative of the target function. */
private transient UnivariateRealFunction derivative;
@ -52,12 +52,12 @@ public class NewtonSolver extends UnivariateRealSolverImpl {
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function or derivative
* @throws IllegalArgumentException if min is not less than max
*/
public double solve(double min, double max) throws ConvergenceException,
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
return solve(min, max, UnivariateRealSolverUtils.midpoint(min, max));
}
@ -69,13 +69,13 @@ public class NewtonSolver extends UnivariateRealSolverImpl {
* @param max the upper bound for the interval (ignored).
* @param startValue the start value to use.
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function or derivative
* @throws IllegalArgumentException if startValue is not between min and max
*/
public double solve(double min, double max, double startValue)
throws ConvergenceException, FunctionEvaluationException {
throws MaxIterationsExceededException, FunctionEvaluationException {
clearResult();
verifySequence(min, startValue, max);
@ -96,8 +96,7 @@ public class NewtonSolver extends UnivariateRealSolverImpl {
++i;
}
throw new ConvergenceException
("Maximum number of iterations exceeded " + i);
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**

316
src/java/org/apache/commons/math/analysis/RiddersSolver.java Executable file → Normal file
View File

@ -1,159 +1,157 @@
/*
* 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;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
* Ridders' Method</a> for root finding of real univariate functions. For
* reference, see C. Ridders, <i>A new algorithm for computing a single root
* of a real continuous function </i>, IEEE Transactions on Circuits and
* Systems, 26 (1979), 979 - 980.
* <p>
* The function should be continuous but not necessarily smooth.
*
* @version $Revision$ $Date$
*/
public class RiddersSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
static final long serialVersionUID = -4703139035737911735L;
/**
* Construct a solver for the given function.
*
* @param f function to solve
*/
public RiddersSolver(UnivariateRealFunction f) {
super(f, 100, 1E-6);
}
/**
* Find a root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
ConvergenceException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a root in the given interval.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws ConvergenceException,
FunctionEvaluationException {
// [x1, x2] is the bracketing interval in each iteration
// x3 is the midpoint of [x1, x2]
// x is the new root approximation and an endpoint of the new interval
double x1, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance;
x1 = min; y1 = f.value(x1);
x2 = max; y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y1 == 0.0) { return min; }
if (y2 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// calculate the new root approximation
x3 = 0.5 * (x1 + x2);
y3 = f.value(x3);
if (Math.abs(y3) <= functionValueAccuracy) {
setResult(x3, i);
return result;
}
delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
(x3 - x1) / Math.sqrt(delta);
x = x3 - correction; // correction != 0
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the new interval for next iteration
// Ridders' method guarantees x1 < x < x2
if (correction > 0.0) { // x1 < x < x3
if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
x2 = x; y2 = y;
} else {
x1 = x; x2 = x3;
y1 = y; y2 = y3;
}
} else { // x3 < x < x2
if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
x1 = x; y1 = y;
} else {
x1 = x3; x2 = x;
y1 = y3; y2 = y;
}
}
oldx = x;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
}
}
/*
* 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;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
* Ridders' Method</a> for root finding of real univariate functions. For
* reference, see C. Ridders, <i>A new algorithm for computing a single root
* of a real continuous function </i>, IEEE Transactions on Circuits and
* Systems, 26 (1979), 979 - 980.
* <p>
* The function should be continuous but not necessarily smooth.
*
* @version $Revision$ $Date$
*/
public class RiddersSolver extends UnivariateRealSolverImpl {
/** serializable version identifier */
private static final long serialVersionUID = -4703139035737911735L;
/**
* Construct a solver for the given function.
*
* @param f function to solve
*/
public RiddersSolver(UnivariateRealFunction f) {
super(f, 100, 1E-6);
}
/**
* Find a root in the given interval with initial value.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max, double initial) throws
MaxIterationsExceededException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(min, initial);
} else {
return solve(initial, max);
}
}
/**
* Find a root in the given interval.
* <p>
* Requires bracketing condition.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
// [x1, x2] is the bracketing interval in each iteration
// x3 is the midpoint of [x1, x2]
// x is the new root approximation and an endpoint of the new interval
double x1, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance;
x1 = min; y1 = f.value(x1);
x2 = max; y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y1 == 0.0) { return min; }
if (y2 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// calculate the new root approximation
x3 = 0.5 * (x1 + x2);
y3 = f.value(x3);
if (Math.abs(y3) <= functionValueAccuracy) {
setResult(x3, i);
return result;
}
delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
(x3 - x1) / Math.sqrt(delta);
x = x3 - correction; // correction != 0
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the new interval for next iteration
// Ridders' method guarantees x1 < x < x2
if (correction > 0.0) { // x1 < x < x3
if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
x2 = x; y2 = y;
} else {
x1 = x; x2 = x3;
y1 = y; y2 = y3;
}
} else { // x3 < x < x2
if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
x1 = x; y1 = y;
} else {
x1 = x3; x2 = x;
y1 = y3; y2 = y;
}
}
oldx = x;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

10
src/java/org/apache/commons/math/analysis/RombergIntegrator.java
View File

@ -16,8 +16,8 @@
*/
package org.apache.commons.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
@ -34,7 +34,7 @@ import org.apache.commons.math.FunctionEvaluationException;
public class RombergIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = -1058849527738180243L;
private static final long serialVersionUID = -1058849527738180243L;
/**
* Construct an integrator for the given function.
@ -51,13 +51,13 @@ public class RombergIntegrator extends UnivariateRealIntegratorImpl {
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
public double integrate(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1, j, m = maximalIterationCount + 1;
@ -89,7 +89,7 @@ public class RombergIntegrator extends UnivariateRealIntegratorImpl {
olds = s;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**

14
src/java/org/apache/commons/math/analysis/SecantSolver.java
View File

@ -18,8 +18,8 @@ package org.apache.commons.math.analysis;
import java.io.Serializable;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
@ -42,7 +42,7 @@ public class SecantSolver extends UnivariateRealSolverImpl implements Serializab
/** Serializable version identifier */
private static final long serialVersionUID = 1984971194738974867L;
/**
* Construct a solver for the given function.
* @param f function to solve.
@ -58,14 +58,14 @@ public class SecantSolver extends UnivariateRealSolverImpl implements Serializab
* @param max the upper bound for the interval
* @param initial the start value to use (ignored)
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
*/
public double solve(double min, double max, double initial)
throws ConvergenceException, FunctionEvaluationException {
throws MaxIterationsExceededException, FunctionEvaluationException {
return solve(min, max);
}
@ -75,13 +75,13 @@ public class SecantSolver extends UnivariateRealSolverImpl implements Serializab
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return the value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
*/
public double solve(double min, double max) throws ConvergenceException,
public double solve(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException {
clearResult();
@ -151,7 +151,7 @@ public class SecantSolver extends UnivariateRealSolverImpl implements Serializab
oldDelta = x2 - x1;
i++;
}
throw new ConvergenceException("Maximal iteration number exceeded" + i);
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

10
src/java/org/apache/commons/math/analysis/SimpsonIntegrator.java
View File

@ -16,8 +16,8 @@
*/
package org.apache.commons.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements the <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
@ -33,7 +33,7 @@ import org.apache.commons.math.FunctionEvaluationException;
public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = 3405465123320678216L;
private static final long serialVersionUID = 3405465123320678216L;
/**
* Construct an integrator for the given function.
@ -50,13 +50,13 @@ public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
public double integrate(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1;
@ -88,7 +88,7 @@ public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
oldt = t;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**

10
src/java/org/apache/commons/math/analysis/TrapezoidIntegrator.java
View File

@ -16,8 +16,8 @@
*/
package org.apache.commons.math.analysis;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html">
@ -32,7 +32,7 @@ import org.apache.commons.math.FunctionEvaluationException;
public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
/** serializable version identifier */
static final long serialVersionUID = 4978222553983172543L;
private static final long serialVersionUID = 4978222553983172543L;
/** intermediate result */
private double s;
@ -91,13 +91,13 @@ public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value of integral
* @throws ConvergenceException if the maximum iteration count is exceeded
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the integrator detects convergence problems otherwise
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double integrate(double min, double max) throws ConvergenceException,
public double integrate(double min, double max) throws MaxIterationsExceededException,
FunctionEvaluationException, IllegalArgumentException {
int i = 1;
@ -119,7 +119,7 @@ public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
oldt = t;
i++;
}
throw new ConvergenceException("Maximum number of iterations exceeded.");
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**

11
src/java/org/apache/commons/math/special/Gamma.java
View File

@ -18,8 +18,8 @@ package org.apache.commons.math.special;
import java.io.Serializable;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.util.ContinuedFraction;
/**
@ -30,6 +30,9 @@ import org.apache.commons.math.util.ContinuedFraction;
*/
public class Gamma implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = -6587513359895466954L;
/** Maximum allowed numerical error. */
private static final double DEFAULT_EPSILON = 10e-9;
@ -174,8 +177,7 @@ public class Gamma implements Serializable {
sum = sum + an;
}
if (n >= maxIterations) {
throw new ConvergenceException(
"maximum number of iterations reached");
throw new MaxIterationsExceededException(maxIterations);
} else {
ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
}
@ -239,6 +241,9 @@ public class Gamma implements Serializable {
} else {
// create continued fraction
ContinuedFraction cf = new ContinuedFraction() {
private static final long serialVersionUID = 5378525034886164398L;
protected double getA(int n, double x) {
return ((2.0 * n) + 1.0) - a + x;
}

10
src/java/org/apache/commons/math/util/ContinuedFraction.java
View File

@ -20,6 +20,7 @@ import java.io.Serializable;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MaxIterationsExceededException;
/**
* Provides a generic means to evaluate continued fractions. Subclasses simply
@ -153,8 +154,8 @@ public abstract class ContinuedFraction implements Serializable {
} else {
// can not scale an convergent is unbounded.
throw new ConvergenceException(
"Continued fraction convergents diverged to +/- " +
"infinity.");
"Continued fraction convergents diverged to +/- infinity for value {0}",
new Object[] { new Double(x) });
}
}
double r = p2 / q2;
@ -169,8 +170,9 @@ public abstract class ContinuedFraction implements Serializable {
}
if (n >= maxIterations) {
throw new ConvergenceException(
"Continued fraction convergents failed to converge.");
throw new MaxIterationsExceededException(maxIterations,
"Continued fraction convergents failed to converge for value {0}",
new Object[] { new Double(x) });
}
return c;

38
src/test/org/apache/commons/math/MaxIterationsExceededExceptionTest.java Normal file
View File

@ -0,0 +1,38 @@
/*
* 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;
import java.util.Locale;
import junit.framework.TestCase;
/**
* @version $Revision:$
*/
public class MaxIterationsExceededExceptionTest extends TestCase {
public void testConstructor(){
MaxIterationsExceededException ex = new MaxIterationsExceededException(1000000);
assertNull(ex.getCause());
assertNotNull(ex.getMessage());
assertTrue(ex.getMessage().indexOf("1,000,000") > 0);
assertEquals(1000000, ex.getMaxIterations());
assertFalse(ex.getMessage().equals(ex.getMessage(Locale.FRENCH)));
}
}