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