src/java/org/apache/commons/math/analysis/LaguerreSolver.java - commons-math - Git at Google

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

 /**
  * 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$ $Date$
  * @since 1.2
  */
 public class LaguerreSolver extends UnivariateRealSolverImpl {

     /** serializable version identifier */
     private static final long serialVersionUID = -3775334783473775723L;

     /** polynomial function to solve */
     private PolynomialFunction p;

     /**
      * Construct a solver for the given function.
      *
      * @param f function to solve
      * @throws IllegalArgumentException if function is not polynomial
      */
     public LaguerreSolver(UnivariateRealFunction f) throws
         IllegalArgumentException {

         super(f, 100, 1E-6);
         if (f instanceof PolynomialFunction) {
             p = (PolynomialFunction)f;
         } else {
             throw new IllegalArgumentException("Function is not polynomial.");
         }
     }

     /**
      * Returns a copy of the polynomial function.
      *
      * @return a fresh copy of the polynomial function
      */
     public PolynomialFunction getPolynomialFunction() {
         return new PolynomialFunction(p.getCoefficients());
     }

     /**
      * Find a real root in the given interval with initial value.
      * <p>
      * Requires bracketing condition.</p>
      *
      * @param min the lower bound for the interval
      * @param max the upper bound for the interval
      * @param initial the start value to use
      * @return the point at which the function value is zero
      * @throws ConvergenceException if the maximum iteration count is exceeded
      * or the solver detects convergence problems otherwise
      * @throws FunctionEvaluationException if an error occurs evaluating the
      * function
      * @throws IllegalArgumentException if any parameters are invalid
      */
     public double solve(double min, double max, double initial) throws
         ConvergenceException, FunctionEvaluationException {

         // check for zeros before verifying bracketing
         if (p.value(min) == 0.0) { return min; }
         if (p.value(max) == 0.0) { return max; }
         if (p.value(initial) == 0.0) { return initial; }

         verifyBracketing(min, max, p);
         verifySequence(min, initial, max);
         if (isBracketing(min, initial, p)) {
             return solve(min, initial);
         } else {
             return solve(initial, max);
         }
     }

     /**
      * Find a real root in the given interval.
      * <p>
      * Despite the bracketing condition, the root returned by solve(Complex[],
      * Complex) may not be a real zero inside [min, max]. For example,
      * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
      * another initial value, or, as we did here, call solveAll() to obtain
      * all roots and pick up the one that we're looking for.</p>
      *
      * @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 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);
     }
 }
	/*
	* 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;

	/**
	* 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$ $Date$
	* @since 1.2
	*/
	public class LaguerreSolver extends UnivariateRealSolverImpl {

	/** serializable version identifier */
	private static final long serialVersionUID = -3775334783473775723L;

	/** polynomial function to solve */
	private PolynomialFunction p;

	/**
	* Construct a solver for the given function.
	*
	* @param f function to solve
	* @throws IllegalArgumentException if function is not polynomial
	*/
	public LaguerreSolver(UnivariateRealFunction f) throws
	IllegalArgumentException {

	super(f, 100, 1E-6);
	if (f instanceof PolynomialFunction) {
	p = (PolynomialFunction)f;
	} else {
	throw new IllegalArgumentException("Function is not polynomial.");
	}
	}

	/**
	* Returns a copy of the polynomial function.
	*
	* @return a fresh copy of the polynomial function
	*/
	public PolynomialFunction getPolynomialFunction() {
	return new PolynomialFunction(p.getCoefficients());
	}

	/**
	* Find a real root in the given interval with initial value.
	* <p>
	* Requires bracketing condition.</p>
	*
	* @param min the lower bound for the interval
	* @param max the upper bound for the interval
	* @param initial the start value to use
	* @return the point at which the function value is zero
	* @throws ConvergenceException if the maximum iteration count is exceeded
	* or the solver detects convergence problems otherwise
	* @throws FunctionEvaluationException if an error occurs evaluating the
	* function
	* @throws IllegalArgumentException if any parameters are invalid
	*/
	public double solve(double min, double max, double initial) throws
	ConvergenceException, FunctionEvaluationException {

	// check for zeros before verifying bracketing
	if (p.value(min) == 0.0) { return min; }
	if (p.value(max) == 0.0) { return max; }
	if (p.value(initial) == 0.0) { return initial; }

	verifyBracketing(min, max, p);
	verifySequence(min, initial, max);
	if (isBracketing(min, initial, p)) {
	return solve(min, initial);
	} else {
	return solve(initial, max);
	}
	}

	/**
	* Find a real root in the given interval.
	* <p>
	* Despite the bracketing condition, the root returned by solve(Complex[],
	* Complex) may not be a real zero inside [min, max]. For example,
	* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
	* another initial value, or, as we did here, call solveAll() to obtain
	* all roots and pick up the one that we're looking for.</p>
	*
	* @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 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);
	}
	}