1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.analysis.solvers;
18
19 import org.apache.commons.math.ConvergenceException;
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.MathRuntimeException;
22 import org.apache.commons.math.MaxIterationsExceededException;
23 import org.apache.commons.math.analysis.UnivariateRealFunction;
24 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
25 import org.apache.commons.math.complex.Complex;
26
27
28
29
30
31
32
33
34
35
36
37
38
39 public class LaguerreSolver extends UnivariateRealSolverImpl {
40
41
42
43 @Deprecated
44 private PolynomialFunction p;
45
46
47
48
49
50
51
52
53
54
55
56 @Deprecated
57 public LaguerreSolver(UnivariateRealFunction f) throws
58 IllegalArgumentException {
59 super(f, 100, 1E-6);
60 if (f instanceof PolynomialFunction) {
61 p = (PolynomialFunction) f;
62 } else {
63 throw MathRuntimeException.createIllegalArgumentException("function is not polynomial");
64 }
65 }
66
67
68
69
70 public LaguerreSolver() {
71 super(100, 1E-6);
72 }
73
74
75
76
77
78
79
80 @Deprecated
81 public PolynomialFunction getPolynomialFunction() {
82 return new PolynomialFunction(p.getCoefficients());
83 }
84
85
86 @Deprecated
87 public double solve(final double min, final double max)
88 throws ConvergenceException, FunctionEvaluationException {
89 return solve(p, min, max);
90 }
91
92
93 @Deprecated
94 public double solve(final double min, final double max, final double initial)
95 throws ConvergenceException, FunctionEvaluationException {
96 return solve(p, min, max, initial);
97 }
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115 public double solve(final UnivariateRealFunction f,
116 final double min, final double max, final double initial)
117 throws ConvergenceException, FunctionEvaluationException {
118
119
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153 public double solve(final UnivariateRealFunction f,
154 final double min, final double max)
155 throws ConvergenceException, FunctionEvaluationException {
156
157
158 if (!(f instanceof PolynomialFunction)) {
159 throw MathRuntimeException.createIllegalArgumentException("function is not polynomial");
160 }
161
162
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
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
189 throw new ConvergenceException();
190 }
191
192
193
194
195
196
197
198
199
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
210
211
212
213
214
215
216
217
218
219
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
234
235
236
237
238
239
240
241
242
243
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];
255 for (int i = 0; i <= n; i++) {
256 c[i] = coefficients[i];
257 }
258
259
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
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
283
284
285
286
287
288
289
290
291
292
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
318
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
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
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
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
354
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 }