Commons Math example source code file (LaguerreSolver.java)

This example Commons Math source code file (LaguerreSolver.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" ^TM.

Java - Commons Math tags/keywords

complex, complex, convergenceexception, convergenceexception, deprecated, functionevaluationexception, functionevaluationexception, g, laguerresolver, maxiterationsexceededexception, non_positive_degree_message, polynomialfunction, polynomialfunction, univariaterealfunction

The Commons Math LaguerreSolver.java source code

/*
 * 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;

/**
 * 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,
 * 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: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $
 * @since 1.2
 */
public class LaguerreSolver extends UnivariateRealSolverImpl {

    /** Message for non-polynomial function. */
    private static final String NON_POLYNOMIAL_FUNCTION_MESSAGE =
        "function is not polynomial";

    /** Message for non-positive degree. */
    private static final String NON_POSITIVE_DEGREE_MESSAGE =
        "polynomial degree must be positive: degree={0}";

    /** 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(NON_POLYNOMIAL_FUNCTION_MESSAGE);
        }
    }

    /**
     * Construct a solver.
     */
    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
     * @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(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
     * @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(final UnivariateRealFunction f,
                        final double min, final double max)
        throws ConvergenceException, FunctionEvaluationException {

        // check function type
        if (!(f instanceof PolynomialFunction)) {
            throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE);
        }

        // 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 = 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 MathRuntimeException.createIllegalArgumentException(
                  NON_POSITIVE_DEGREE_MESSAGE, 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 MathRuntimeException.createIllegalArgumentException(
                  NON_POSITIVE_DEGREE_MESSAGE, 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 = 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 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);
    }
}