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.ConvergenceException;
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.MathRuntimeException;
022 import org.apache.commons.math.MaxIterationsExceededException;
023 import org.apache.commons.math.analysis.UnivariateRealFunction;
024 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
025 import org.apache.commons.math.complex.Complex;
026
027 /**
028 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
029 * Laguerre's Method</a> for root finding of real coefficient polynomials.
030 * For reference, see <b>A First Course in Numerical Analysis</b>,
031 * ISBN 048641454X, chapter 8.
032 * <p>
033 * Laguerre's method is global in the sense that it can start with any initial
034 * approximation and be able to solve all roots from that point.</p>
035 *
036 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
037 * @since 1.2
038 */
039 public class LaguerreSolver extends UnivariateRealSolverImpl {
040 /** polynomial function to solve.
041 * @deprecated as of 2.0 the function is not stored anymore in the instance
042 */
043 @Deprecated
044 private PolynomialFunction p;
045
046 /**
047 * Construct a solver for the given function.
048 *
049 * @param f function to solve
050 * @throws IllegalArgumentException if function is not polynomial
051 * @deprecated as of 2.0 the function to solve is passed as an argument
052 * to the {@link #solve(UnivariateRealFunction, double, double)} or
053 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
054 * method.
055 */
056 @Deprecated
057 public LaguerreSolver(UnivariateRealFunction f) throws
058 IllegalArgumentException {
059 super(f, 100, 1E-6);
060 if (f instanceof PolynomialFunction) {
061 p = (PolynomialFunction) f;
062 } else {
063 throw MathRuntimeException.createIllegalArgumentException("function is not polynomial");
064 }
065 }
066
067 /**
068 * Construct a solver.
069 */
070 public LaguerreSolver() {
071 super(100, 1E-6);
072 }
073
074 /**
075 * Returns a copy of the polynomial function.
076 *
077 * @return a fresh copy of the polynomial function
078 * @deprecated as of 2.0 the function is not stored anymore within the instance.
079 */
080 @Deprecated
081 public PolynomialFunction getPolynomialFunction() {
082 return new PolynomialFunction(p.getCoefficients());
083 }
084
085 /** {@inheritDoc} */
086 @Deprecated
087 public double solve(final double min, final double max)
088 throws ConvergenceException, FunctionEvaluationException {
089 return solve(p, min, max);
090 }
091
092 /** {@inheritDoc} */
093 @Deprecated
094 public double solve(final double min, final double max, final double initial)
095 throws ConvergenceException, FunctionEvaluationException {
096 return solve(p, min, max, initial);
097 }
098
099 /**
100 * Find a real root in the given interval with initial value.
101 * <p>
102 * Requires bracketing condition.</p>
103 *
104 * @param f function to solve (must be polynomial)
105 * @param min the lower bound for the interval
106 * @param max the upper bound for the interval
107 * @param initial the start value to use
108 * @return the point at which the function value is zero
109 * @throws ConvergenceException if the maximum iteration count is exceeded
110 * or the solver detects convergence problems otherwise
111 * @throws FunctionEvaluationException if an error occurs evaluating the
112 * function
113 * @throws IllegalArgumentException if any parameters are invalid
114 */
115 public double solve(final UnivariateRealFunction f,
116 final double min, final double max, final double initial)
117 throws ConvergenceException, FunctionEvaluationException {
118
119 // check for zeros before verifying bracketing
120 if (f.value(min) == 0.0) { return min; }
121 if (f.value(max) == 0.0) { return max; }
122 if (f.value(initial) == 0.0) { return initial; }
123
124 verifyBracketing(min, max, f);
125 verifySequence(min, initial, max);
126 if (isBracketing(min, initial, f)) {
127 return solve(f, min, initial);
128 } else {
129 return solve(f, initial, max);
130 }
131
132 }
133
134 /**
135 * Find a real root in the given interval.
136 * <p>
137 * Despite the bracketing condition, the root returned by solve(Complex[],
138 * Complex) may not be a real zero inside [min, max]. For example,
139 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
140 * another initial value, or, as we did here, call solveAll() to obtain
141 * all roots and pick up the one that we're looking for.</p>
142 *
143 * @param f the function to solve
144 * @param min the lower bound for the interval
145 * @param max the upper bound for the interval
146 * @return the point at which the function value is zero
147 * @throws ConvergenceException if the maximum iteration count is exceeded
148 * or the solver detects convergence problems otherwise
149 * @throws FunctionEvaluationException if an error occurs evaluating the
150 * function
151 * @throws IllegalArgumentException if any parameters are invalid
152 */
153 public double solve(final UnivariateRealFunction f,
154 final double min, final double max)
155 throws ConvergenceException, FunctionEvaluationException {
156
157 // check function type
158 if (!(f instanceof PolynomialFunction)) {
159 throw MathRuntimeException.createIllegalArgumentException("function is not polynomial");
160 }
161
162 // check for zeros before verifying bracketing
163 if (f.value(min) == 0.0) { return min; }
164 if (f.value(max) == 0.0) { return max; }
165 verifyBracketing(min, max, f);
166
167 double coefficients[] = ((PolynomialFunction) f).getCoefficients();
168 Complex c[] = new Complex[coefficients.length];
169 for (int i = 0; i < coefficients.length; i++) {
170 c[i] = new Complex(coefficients[i], 0.0);
171 }
172 Complex initial = new Complex(0.5 * (min + max), 0.0);
173 Complex z = solve(c, initial);
174 if (isRootOK(min, max, z)) {
175 setResult(z.getReal(), iterationCount);
176 return result;
177 }
178
179 // solve all roots and select the one we're seeking
180 Complex[] root = solveAll(c, initial);
181 for (int i = 0; i < root.length; i++) {
182 if (isRootOK(min, max, root[i])) {
183 setResult(root[i].getReal(), iterationCount);
184 return result;
185 }
186 }
187
188 // should never happen
189 throw new ConvergenceException();
190 }
191
192 /**
193 * Returns true iff the given complex root is actually a real zero
194 * in the given interval, within the solver tolerance level.
195 *
196 * @param min the lower bound for the interval
197 * @param max the upper bound for the interval
198 * @param z the complex root
199 * @return true iff z is the sought-after real zero
200 */
201 protected boolean isRootOK(double min, double max, Complex z) {
202 double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);
203 return (isSequence(min, z.getReal(), max)) &&
204 (Math.abs(z.getImaginary()) <= tolerance ||
205 z.abs() <= functionValueAccuracy);
206 }
207
208 /**
209 * Find all complex roots for the polynomial with the given coefficients,
210 * starting from the given initial value.
211 *
212 * @param coefficients the polynomial coefficients array
213 * @param initial the start value to use
214 * @return the point at which the function value is zero
215 * @throws ConvergenceException if the maximum iteration count is exceeded
216 * or the solver detects convergence problems otherwise
217 * @throws FunctionEvaluationException if an error occurs evaluating the
218 * function
219 * @throws IllegalArgumentException if any parameters are invalid
220 */
221 public Complex[] solveAll(double coefficients[], double initial) throws
222 ConvergenceException, FunctionEvaluationException {
223
224 Complex c[] = new Complex[coefficients.length];
225 Complex z = new Complex(initial, 0.0);
226 for (int i = 0; i < c.length; i++) {
227 c[i] = new Complex(coefficients[i], 0.0);
228 }
229 return solveAll(c, z);
230 }
231
232 /**
233 * Find all complex roots for the polynomial with the given coefficients,
234 * starting from the given initial value.
235 *
236 * @param coefficients the polynomial coefficients array
237 * @param initial the start value to use
238 * @return the point at which the function value is zero
239 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
240 * or the solver detects convergence problems otherwise
241 * @throws FunctionEvaluationException if an error occurs evaluating the
242 * function
243 * @throws IllegalArgumentException if any parameters are invalid
244 */
245 public Complex[] solveAll(Complex coefficients[], Complex initial) throws
246 MaxIterationsExceededException, FunctionEvaluationException {
247
248 int n = coefficients.length - 1;
249 int iterationCount = 0;
250 if (n < 1) {
251 throw MathRuntimeException.createIllegalArgumentException(
252 "polynomial degree must be positive: degree={0}", n);
253 }
254 Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial
255 for (int i = 0; i <= n; i++) {
256 c[i] = coefficients[i];
257 }
258
259 // solve individual root successively
260 Complex root[] = new Complex[n];
261 for (int i = 0; i < n; i++) {
262 Complex subarray[] = new Complex[n-i+1];
263 System.arraycopy(c, 0, subarray, 0, subarray.length);
264 root[i] = solve(subarray, initial);
265 // polynomial deflation using synthetic division
266 Complex newc = c[n-i];
267 Complex oldc = null;
268 for (int j = n-i-1; j >= 0; j--) {
269 oldc = c[j];
270 c[j] = newc;
271 newc = oldc.add(newc.multiply(root[i]));
272 }
273 iterationCount += this.iterationCount;
274 }
275
276 resultComputed = true;
277 this.iterationCount = iterationCount;
278 return root;
279 }
280
281 /**
282 * Find a complex root for the polynomial with the given coefficients,
283 * starting from the given initial value.
284 *
285 * @param coefficients the polynomial coefficients array
286 * @param initial the start value to use
287 * @return the point at which the function value is zero
288 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
289 * or the solver detects convergence problems otherwise
290 * @throws FunctionEvaluationException if an error occurs evaluating the
291 * function
292 * @throws IllegalArgumentException if any parameters are invalid
293 */
294 public Complex solve(Complex coefficients[], Complex initial) throws
295 MaxIterationsExceededException, FunctionEvaluationException {
296
297 int n = coefficients.length - 1;
298 if (n < 1) {
299 throw MathRuntimeException.createIllegalArgumentException(
300 "polynomial degree must be positive: degree={0}", n);
301 }
302 Complex N = new Complex(n, 0.0);
303 Complex N1 = new Complex((n-1), 0.0);
304
305 int i = 1;
306 Complex pv = null;
307 Complex dv = null;
308 Complex d2v = null;
309 Complex G = null;
310 Complex G2 = null;
311 Complex H = null;
312 Complex delta = null;
313 Complex denominator = null;
314 Complex z = initial;
315 Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
316 while (i <= maximalIterationCount) {
317 // Compute pv (polynomial value), dv (derivative value), and
318 // d2v (second derivative value) simultaneously.
319 pv = coefficients[n];
320 dv = Complex.ZERO;
321 d2v = Complex.ZERO;
322 for (int j = n-1; j >= 0; j--) {
323 d2v = dv.add(z.multiply(d2v));
324 dv = pv.add(z.multiply(dv));
325 pv = coefficients[j].add(z.multiply(pv));
326 }
327 d2v = d2v.multiply(new Complex(2.0, 0.0));
328
329 // check for convergence
330 double tolerance = Math.max(relativeAccuracy * z.abs(),
331 absoluteAccuracy);
332 if ((z.subtract(oldz)).abs() <= tolerance) {
333 resultComputed = true;
334 iterationCount = i;
335 return z;
336 }
337 if (pv.abs() <= functionValueAccuracy) {
338 resultComputed = true;
339 iterationCount = i;
340 return z;
341 }
342
343 // now pv != 0, calculate the new approximation
344 G = dv.divide(pv);
345 G2 = G.multiply(G);
346 H = G2.subtract(d2v.divide(pv));
347 delta = N1.multiply((N.multiply(H)).subtract(G2));
348 // choose a denominator larger in magnitude
349 Complex deltaSqrt = delta.sqrt();
350 Complex dplus = G.add(deltaSqrt);
351 Complex dminus = G.subtract(deltaSqrt);
352 denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
353 // Perturb z if denominator is zero, for instance,
354 // p(x) = x^3 + 1, z = 0.
355 if (denominator.equals(new Complex(0.0, 0.0))) {
356 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
357 oldz = new Complex(Double.POSITIVE_INFINITY,
358 Double.POSITIVE_INFINITY);
359 } else {
360 oldz = z;
361 z = z.subtract(N.divide(denominator));
362 }
363 i++;
364 }
365 throw new MaxIterationsExceededException(maximalIterationCount);
366 }
367 }