- * 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. - *
- * 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. - *
- * 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 + * Laguerre's Method for root finding of real coefficient polynomials. + * For reference, see A First Course in Numerical Analysis, + * ISBN 048641454X, chapter 8. + *
+ * 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. + *
+ * 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. + *
+ * 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); + } +} diff --git a/src/java/org/apache/commons/math/analysis/MullerSolver.java b/src/java/org/apache/commons/math/analysis/MullerSolver.java old mode 100755 new mode 100644 index dad3a5ab0..5d1e765e4 --- a/src/java/org/apache/commons/math/analysis/MullerSolver.java +++ b/src/java/org/apache/commons/math/analysis/MullerSolver.java @@ -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 - * Muller's Method for root finding of real univariate functions. For - * reference, see Elementary Numerical Analysis, ISBN 0070124477, - * chapter 3. - *
- * 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. - *
- * 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. - *
- * 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. - *
- * 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. - *
- * 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. - *
- * 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. - *
- * 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. - *
- * 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 + * Muller's Method for root finding of real univariate functions. For + * reference, see Elementary Numerical Analysis, ISBN 0070124477, + * chapter 3. + *
+ * 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. + *
+ * 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. + *
+ * 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. + *
+ * 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. + *
+ * 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. + *
+ * 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. + *
+ * 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. + *
+ * 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); + } +} diff --git a/src/java/org/apache/commons/math/analysis/NewtonSolver.java b/src/java/org/apache/commons/math/analysis/NewtonSolver.java index 934e817c6..ea9f5a7fd 100644 --- a/src/java/org/apache/commons/math/analysis/NewtonSolver.java +++ b/src/java/org/apache/commons/math/analysis/NewtonSolver.java @@ -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 @@ -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); } /** diff --git a/src/java/org/apache/commons/math/analysis/RiddersSolver.java b/src/java/org/apache/commons/math/analysis/RiddersSolver.java old mode 100755 new mode 100644 index b12b02afe..298597cfc --- a/src/java/org/apache/commons/math/analysis/RiddersSolver.java +++ b/src/java/org/apache/commons/math/analysis/RiddersSolver.java @@ -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 - * Ridders' Method for root finding of real univariate functions. For - * reference, see C. Ridders, A new algorithm for computing a single root - * of a real continuous function , IEEE Transactions on Circuits and - * Systems, 26 (1979), 979 - 980. - *
- * 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. - *
- * 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. - *
- * 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 + * Ridders' Method for root finding of real univariate functions. For + * reference, see C. Ridders, A new algorithm for computing a single root + * of a real continuous function , IEEE Transactions on Circuits and + * Systems, 26 (1979), 979 - 980. + *
+ * 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. + *
+ * 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. + *
+ * 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); + } +} diff --git a/src/java/org/apache/commons/math/analysis/RombergIntegrator.java b/src/java/org/apache/commons/math/analysis/RombergIntegrator.java index 5cf90a019..863f7fd94 100644 --- a/src/java/org/apache/commons/math/analysis/RombergIntegrator.java +++ b/src/java/org/apache/commons/math/analysis/RombergIntegrator.java @@ -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 @@ -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); } /** diff --git a/src/java/org/apache/commons/math/analysis/SecantSolver.java b/src/java/org/apache/commons/math/analysis/SecantSolver.java index 1d5fc74f4..862408684 100644 --- a/src/java/org/apache/commons/math/analysis/SecantSolver.java +++ b/src/java/org/apache/commons/math/analysis/SecantSolver.java @@ -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); } } diff --git a/src/java/org/apache/commons/math/analysis/SimpsonIntegrator.java b/src/java/org/apache/commons/math/analysis/SimpsonIntegrator.java index 41d575741..f74dbf6a3 100644 --- a/src/java/org/apache/commons/math/analysis/SimpsonIntegrator.java +++ b/src/java/org/apache/commons/math/analysis/SimpsonIntegrator.java @@ -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 @@ -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); } /** diff --git a/src/java/org/apache/commons/math/analysis/TrapezoidIntegrator.java b/src/java/org/apache/commons/math/analysis/TrapezoidIntegrator.java index 9b2b04398..2f32587a2 100644 --- a/src/java/org/apache/commons/math/analysis/TrapezoidIntegrator.java +++ b/src/java/org/apache/commons/math/analysis/TrapezoidIntegrator.java @@ -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 @@ -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); } /** diff --git a/src/java/org/apache/commons/math/special/Gamma.java b/src/java/org/apache/commons/math/special/Gamma.java index 8c13d205b..8c565cbff 100644 --- a/src/java/org/apache/commons/math/special/Gamma.java +++ b/src/java/org/apache/commons/math/special/Gamma.java @@ -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; } diff --git a/src/java/org/apache/commons/math/util/ContinuedFraction.java b/src/java/org/apache/commons/math/util/ContinuedFraction.java index ef87e3c6b..427d8f16a 100644 --- a/src/java/org/apache/commons/math/util/ContinuedFraction.java +++ b/src/java/org/apache/commons/math/util/ContinuedFraction.java @@ -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; diff --git a/src/test/org/apache/commons/math/MaxIterationsExceededExceptionTest.java b/src/test/org/apache/commons/math/MaxIterationsExceededExceptionTest.java new file mode 100644 index 000000000..785640055 --- /dev/null +++ b/src/test/org/apache/commons/math/MaxIterationsExceededExceptionTest.java @@ -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))); + } + +}