001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.analysis.solvers;
018
019 import org.apache.commons.math.complex.Complex;
020 import org.apache.commons.math.exception.NoBracketingException;
021 import org.apache.commons.math.exception.NullArgumentException;
022 import org.apache.commons.math.exception.NoDataException;
023 import org.apache.commons.math.exception.util.LocalizedFormats;
024 import org.apache.commons.math.util.FastMath;
025
026 /**
027 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
028 * Laguerre's Method</a> for root finding of real coefficient polynomials.
029 * For reference, see
030 * <quote>
031 * <b>A First Course in Numerical Analysis</b>
032 * ISBN 048641454X, chapter 8.
033 * </quote>
034 * Laguerre's method is global in the sense that it can start with any initial
035 * approximation and be able to solve all roots from that point.
036 * The algorithm requires a bracketing condition.
037 *
038 * @version $Id: LaguerreSolver.java 1145746 2011-07-12 20:15:01Z psteitz $
039 * @since 1.2
040 */
041 public class LaguerreSolver extends AbstractPolynomialSolver {
042 /** Default absolute accuracy. */
043 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
044 /** Complex solver. */
045 protected ComplexSolver complexSolver = new ComplexSolver();
046
047 /**
048 * Construct a solver with default accuracy (1e-6).
049 */
050 public LaguerreSolver() {
051 this(DEFAULT_ABSOLUTE_ACCURACY);
052 }
053 /**
054 * Construct a solver.
055 *
056 * @param absoluteAccuracy Absolute accuracy.
057 */
058 public LaguerreSolver(double absoluteAccuracy) {
059 super(absoluteAccuracy);
060 }
061 /**
062 * Construct a solver.
063 *
064 * @param relativeAccuracy Relative accuracy.
065 * @param absoluteAccuracy Absolute accuracy.
066 */
067 public LaguerreSolver(double relativeAccuracy,
068 double absoluteAccuracy) {
069 super(relativeAccuracy, absoluteAccuracy);
070 }
071 /**
072 * Construct a solver.
073 *
074 * @param relativeAccuracy Relative accuracy.
075 * @param absoluteAccuracy Absolute accuracy.
076 * @param functionValueAccuracy Function value accuracy.
077 */
078 public LaguerreSolver(double relativeAccuracy,
079 double absoluteAccuracy,
080 double functionValueAccuracy) {
081 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
082 }
083
084 /**
085 * {@inheritDoc}
086 */
087 @Override
088 public double doSolve() {
089 double min = getMin();
090 double max = getMax();
091 double initial = getStartValue();
092 final double functionValueAccuracy = getFunctionValueAccuracy();
093
094 verifySequence(min, initial, max);
095
096 // Return the initial guess if it is good enough.
097 double yInitial = computeObjectiveValue(initial);
098 if (FastMath.abs(yInitial) <= functionValueAccuracy) {
099 return initial;
100 }
101
102 // Return the first endpoint if it is good enough.
103 double yMin = computeObjectiveValue(min);
104 if (FastMath.abs(yMin) <= functionValueAccuracy) {
105 return min;
106 }
107
108 // Reduce interval if min and initial bracket the root.
109 if (yInitial * yMin < 0) {
110 return laguerre(min, initial, yMin, yInitial);
111 }
112
113 // Return the second endpoint if it is good enough.
114 double yMax = computeObjectiveValue(max);
115 if (FastMath.abs(yMax) <= functionValueAccuracy) {
116 return max;
117 }
118
119 // Reduce interval if initial and max bracket the root.
120 if (yInitial * yMax < 0) {
121 return laguerre(initial, max, yInitial, yMax);
122 }
123
124 throw new NoBracketingException(min, max, yMin, yMax);
125 }
126
127 /**
128 * Find a real root in the given interval.
129 *
130 * Despite the bracketing condition, the root returned by
131 * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
132 * not be a real zero inside {@code [min, max]}.
133 * For example, <code>p(x) = x<sup>3</sup> + 1,</code>
134 * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
135 * When it occurs, this code calls
136 * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
137 * in order to obtain all roots and picks up one real root.
138 *
139 * @param lo Lower bound of the search interval.
140 * @param hi Higher bound of the search interval.
141 * @param fLo Function value at the lower bound of the search interval.
142 * @param fHi Function value at the higher bound of the search interval.
143 * @return the point at which the function value is zero.
144 */
145 public double laguerre(double lo, double hi,
146 double fLo, double fHi) {
147 double coefficients[] = getCoefficients();
148 Complex c[] = new Complex[coefficients.length];
149 for (int i = 0; i < coefficients.length; i++) {
150 c[i] = new Complex(coefficients[i], 0);
151 }
152 Complex initial = new Complex(0.5 * (lo + hi), 0);
153 Complex z = complexSolver.solve(c, initial);
154 if (complexSolver.isRoot(lo, hi, z)) {
155 return z.getReal();
156 } else {
157 double r = Double.NaN;
158 // Solve all roots and select the one we are seeking.
159 Complex[] root = complexSolver.solveAll(c, initial);
160 for (int i = 0; i < root.length; i++) {
161 if (complexSolver.isRoot(lo, hi, root[i])) {
162 r = root[i].getReal();
163 break;
164 }
165 }
166 return r;
167 }
168 }
169
170 /**
171 * Class for searching all (complex) roots.
172 */
173 private class ComplexSolver {
174 /**
175 * Check whether the given complex root is actually a real zero
176 * in the given interval, within the solver tolerance level.
177 *
178 * @param min Lower bound for the interval.
179 * @param max Upper bound for the interval.
180 * @param z Complex root.
181 * @return {@code true} if z is a real zero.
182 */
183 public boolean isRoot(double min, double max, Complex z) {
184 if (isSequence(min, z.getReal(), max)) {
185 double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
186 return (FastMath.abs(z.getImaginary()) <= tolerance) ||
187 (z.abs() <= getFunctionValueAccuracy());
188 }
189 return false;
190 }
191
192 /**
193 * Find all complex roots for the polynomial with the given
194 * coefficients, starting from the given initial value.
195 *
196 * @param coefficients Polynomial coefficients.
197 * @param initial Start value.
198 * @return the point at which the function value is zero.
199 * @throws org.apache.commons.math.exception.TooManyEvaluationsException
200 * if the maximum number of evaluations is exceeded.
201 * @throws NullArgumentException if the {@code coefficients} is
202 * {@code null}.
203 * @throws NoDataException if the {@code coefficients} array is empty.
204 */
205 public Complex[] solveAll(Complex coefficients[], Complex initial) {
206 if (coefficients == null) {
207 throw new NullArgumentException();
208 }
209 int n = coefficients.length - 1;
210 if (n == 0) {
211 throw new NoDataException(LocalizedFormats.POLYNOMIAL);
212 }
213 // Coefficients for deflated polynomial.
214 Complex c[] = new Complex[n + 1];
215 for (int i = 0; i <= n; i++) {
216 c[i] = coefficients[i];
217 }
218
219 // Solve individual roots successively.
220 Complex root[] = new Complex[n];
221 for (int i = 0; i < n; i++) {
222 Complex subarray[] = new Complex[n - i + 1];
223 System.arraycopy(c, 0, subarray, 0, subarray.length);
224 root[i] = solve(subarray, initial);
225 // Polynomial deflation using synthetic division.
226 Complex newc = c[n - i];
227 Complex oldc = null;
228 for (int j = n - i - 1; j >= 0; j--) {
229 oldc = c[j];
230 c[j] = newc;
231 newc = oldc.add(newc.multiply(root[i]));
232 }
233 }
234
235 return root;
236 }
237
238 /**
239 * Find a complex root for the polynomial with the given coefficients,
240 * starting from the given initial value.
241 *
242 * @param coefficients Polynomial coefficients.
243 * @param initial Start value.
244 * @return the point at which the function value is zero.
245 * @throws org.apache.commons.math.exception.TooManyEvaluationsException
246 * if the maximum number of evaluations is exceeded.
247 * @throws NullArgumentException if the {@code coefficients} is
248 * {@code null}.
249 * @throws NoDataException if the {@code coefficients} array is empty.
250 */
251 public Complex solve(Complex coefficients[], Complex initial) {
252 if (coefficients == null) {
253 throw new NullArgumentException();
254 }
255
256 int n = coefficients.length - 1;
257 if (n == 0) {
258 throw new NoDataException(LocalizedFormats.POLYNOMIAL);
259 }
260
261 final double absoluteAccuracy = getAbsoluteAccuracy();
262 final double relativeAccuracy = getRelativeAccuracy();
263 final double functionValueAccuracy = getFunctionValueAccuracy();
264
265 Complex N = new Complex(n, 0.0);
266 Complex N1 = new Complex(n - 1, 0.0);
267
268 Complex pv = null;
269 Complex dv = null;
270 Complex d2v = null;
271 Complex G = null;
272 Complex G2 = null;
273 Complex H = null;
274 Complex delta = null;
275 Complex denominator = null;
276 Complex z = initial;
277 Complex oldz = new Complex(Double.POSITIVE_INFINITY,
278 Double.POSITIVE_INFINITY);
279 while (true) {
280 // Compute pv (polynomial value), dv (derivative value), and
281 // d2v (second derivative value) simultaneously.
282 pv = coefficients[n];
283 dv = Complex.ZERO;
284 d2v = Complex.ZERO;
285 for (int j = n-1; j >= 0; j--) {
286 d2v = dv.add(z.multiply(d2v));
287 dv = pv.add(z.multiply(dv));
288 pv = coefficients[j].add(z.multiply(pv));
289 }
290 d2v = d2v.multiply(new Complex(2.0, 0.0));
291
292 // check for convergence
293 double tolerance = FastMath.max(relativeAccuracy * z.abs(),
294 absoluteAccuracy);
295 if ((z.subtract(oldz)).abs() <= tolerance) {
296 return z;
297 }
298 if (pv.abs() <= functionValueAccuracy) {
299 return z;
300 }
301
302 // now pv != 0, calculate the new approximation
303 G = dv.divide(pv);
304 G2 = G.multiply(G);
305 H = G2.subtract(d2v.divide(pv));
306 delta = N1.multiply((N.multiply(H)).subtract(G2));
307 // choose a denominator larger in magnitude
308 Complex deltaSqrt = delta.sqrt();
309 Complex dplus = G.add(deltaSqrt);
310 Complex dminus = G.subtract(deltaSqrt);
311 denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
312 // Perturb z if denominator is zero, for instance,
313 // p(x) = x^3 + 1, z = 0.
314 if (denominator.equals(new Complex(0.0, 0.0))) {
315 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
316 oldz = new Complex(Double.POSITIVE_INFINITY,
317 Double.POSITIVE_INFINITY);
318 } else {
319 oldz = z;
320 z = z.subtract(N.divide(denominator));
321 }
322 incrementEvaluationCount();
323 }
324 }
325 }
326 }