diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java | 434 |
1 files changed, 434 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java b/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java new file mode 100644 index 0000000..e6d8bd8 --- /dev/null +++ b/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java @@ -0,0 +1,434 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.analysis.solvers; + +import org.apache.commons.math.ConvergenceException; +import org.apache.commons.math.FunctionEvaluationException; +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.MaxIterationsExceededException; +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math.complex.Complex; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; + +/** + * 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.</p> + * + * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ + * @since 1.2 + */ +public class LaguerreSolver extends UnivariateRealSolverImpl { + + /** polynomial function to solve. + * @deprecated as of 2.0 the function is not stored anymore in the instance + */ + @Deprecated + private final PolynomialFunction p; + + /** + * Construct a solver for the given function. + * + * @param f function to solve + * @throws IllegalArgumentException if function is not polynomial + * @deprecated as of 2.0 the function to solve is passed as an argument + * to the {@link #solve(UnivariateRealFunction, double, double)} or + * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)} + * method. + */ + @Deprecated + public LaguerreSolver(UnivariateRealFunction f) throws IllegalArgumentException { + super(f, 100, 1E-6); + if (f instanceof PolynomialFunction) { + p = (PolynomialFunction) f; + } else { + throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL); + } + } + + /** + * Construct a solver. + * @deprecated in 2.2 (to be removed in 3.0) + */ + @Deprecated + public LaguerreSolver() { + super(100, 1E-6); + p = null; + } + + /** + * Returns a copy of the polynomial function. + * + * @return a fresh copy of the polynomial function + * @deprecated as of 2.0 the function is not stored anymore within the instance. + */ + @Deprecated + public PolynomialFunction getPolynomialFunction() { + return new PolynomialFunction(p.getCoefficients()); + } + + /** {@inheritDoc} */ + @Deprecated + public double solve(final double min, final double max) + throws ConvergenceException, FunctionEvaluationException { + return solve(p, min, max); + } + + /** {@inheritDoc} */ + @Deprecated + public double solve(final double min, final double max, final double initial) + throws ConvergenceException, FunctionEvaluationException { + return solve(p, min, max, initial); + } + + /** + * Find a real root in the given interval with initial value. + * <p> + * Requires bracketing condition.</p> + * + * @param f function to solve (must be polynomial) + * @param min the lower bound for the interval + * @param max the upper bound for the interval + * @param initial the start value to use + * @param maxEval Maximum number of evaluations. + * @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 + */ + @Override + public double solve(int maxEval, final UnivariateRealFunction f, + final double min, final double max, final double initial) + throws ConvergenceException, FunctionEvaluationException { + setMaximalIterationCount(maxEval); + return solve(f, min, max, initial); + } + + /** + * Find a real root in the given interval with initial value. + * <p> + * Requires bracketing condition.</p> + * + * @param f function to solve (must be polynomial) + * @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 + * @deprecated in 2.2 (to be removed in 3.0). + */ + @Deprecated + public double solve(final UnivariateRealFunction f, + final double min, final double max, final 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(f, min, initial); + } else { + return solve(f, 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.</p> + * + * @param f the function to solve + * @param min the lower bound for the interval + * @param max the upper bound for the interval + * @param maxEval Maximum number of evaluations. + * @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 + */ + @Override + public double solve(int maxEval, final UnivariateRealFunction f, + final double min, final double max) + throws ConvergenceException, FunctionEvaluationException { + setMaximalIterationCount(maxEval); + return solve(f, min, 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.</p> + * + * @param f the function to solve + * @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 + * @deprecated in 2.2 (to be removed in 3.0). + */ + @Deprecated + public double solve(final UnivariateRealFunction f, + final double min, final double max) + throws ConvergenceException, FunctionEvaluationException { + + // check function type + if (!(f instanceof PolynomialFunction)) { + throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL); + } + + // check for zeros before verifying bracketing + if (f.value(min) == 0.0) { return min; } + if (f.value(max) == 0.0) { return max; } + verifyBracketing(min, max, f); + + double coefficients[] = ((PolynomialFunction) f).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 = FastMath.max(relativeAccuracy * z.abs(), absoluteAccuracy); + return (isSequence(min, z.getReal(), max)) && + (FastMath.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 + * @deprecated in 2.2. + */ + @Deprecated + 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 + * @deprecated in 2.2. + */ + @Deprecated + public Complex[] solveAll(Complex coefficients[], Complex initial) throws + MaxIterationsExceededException, FunctionEvaluationException { + + int n = coefficients.length - 1; + int iterationCount = 0; + if (n < 1) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.NON_POSITIVE_POLYNOMIAL_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 + * @deprecated in 2.2. + */ + @Deprecated + public Complex solve(Complex coefficients[], Complex initial) throws + MaxIterationsExceededException, FunctionEvaluationException { + + int n = coefficients.length - 1; + if (n < 1) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n); + } + Complex N = new Complex(n, 0.0); + Complex N1 = new Complex(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 = FastMath.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 deltaSqrt = delta.sqrt(); + Complex dplus = G.add(deltaSqrt); + Complex dminus = G.subtract(deltaSqrt); + 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); + } +} |