1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.analysis.solvers;
18
19 import org.apache.commons.math4.legacy.exception.NoBracketingException;
20 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
21 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
22 import org.apache.commons.math4.core.jdkmath.JdkMath;
23
24 /**
25 * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
26 * Muller's Method</a> for root finding of real univariate functions. For
27 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
28 * chapter 3.
29 * <p>
30 * Muller's method applies to both real and complex functions, but here we
31 * restrict ourselves to real functions.
32 * This class differs from {@link MullerSolver} in the way it avoids complex
33 * operations.</p><p>
34 * Except for the initial [min, max], it does not require bracketing
35 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If a complex
36 * number arises in the computation, we simply use its modulus as a real
37 * approximation.</p>
38 * <p>
39 * Because the interval may not be bracketing, the bisection alternative is
40 * not applicable here. However in practice our treatment usually works
41 * well, especially near real zeroes where the imaginary part of the complex
42 * approximation is often negligible.</p>
43 * <p>
44 * The formulas here do not use divided differences directly.</p>
45 *
46 * @since 1.2
47 * @see MullerSolver
48 */
49 public class MullerSolver2 extends AbstractUnivariateSolver {
50
51 /** Default absolute accuracy. */
52 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
53
54 /**
55 * Construct a solver with default accuracy (1e-6).
56 */
57 public MullerSolver2() {
58 this(DEFAULT_ABSOLUTE_ACCURACY);
59 }
60 /**
61 * Construct a solver.
62 *
63 * @param absoluteAccuracy Absolute accuracy.
64 */
65 public MullerSolver2(double absoluteAccuracy) {
66 super(absoluteAccuracy);
67 }
68 /**
69 * Construct a solver.
70 *
71 * @param relativeAccuracy Relative accuracy.
72 * @param absoluteAccuracy Absolute accuracy.
73 */
74 public MullerSolver2(double relativeAccuracy,
75 double absoluteAccuracy) {
76 super(relativeAccuracy, absoluteAccuracy);
77 }
78
79 /**
80 * {@inheritDoc}
81 */
82 @Override
83 protected double doSolve()
84 throws TooManyEvaluationsException,
85 NumberIsTooLargeException,
86 NoBracketingException {
87 final double min = getMin();
88 final double max = getMax();
89
90 verifyInterval(min, max);
91
92 final double relativeAccuracy = getRelativeAccuracy();
93 final double absoluteAccuracy = getAbsoluteAccuracy();
94 final double functionValueAccuracy = getFunctionValueAccuracy();
95
96 // x2 is the last root approximation
97 // x is the new approximation and new x2 for next round
98 // x0 < x1 < x2 does not hold here
99
100 double x0 = min;
101 double y0 = computeObjectiveValue(x0);
102 if (JdkMath.abs(y0) < functionValueAccuracy) {
103 return x0;
104 }
105 double x1 = max;
106 double y1 = computeObjectiveValue(x1);
107 if (JdkMath.abs(y1) < functionValueAccuracy) {
108 return x1;
109 }
110
111 if(y0 * y1 > 0) {
112 throw new NoBracketingException(x0, x1, y0, y1);
113 }
114
115 double x2 = 0.5 * (x0 + x1);
116 double y2 = computeObjectiveValue(x2);
117
118 double oldx = Double.POSITIVE_INFINITY;
119 while (true) {
120 // quadratic interpolation through x0, x1, x2
121 final double q = (x2 - x1) / (x1 - x0);
122 final double a = q * (y2 - (1 + q) * y1 + q * y0);
123 final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
124 final double c = (1 + q) * y2;
125 final double delta = b * b - 4 * a * c;
126 double x;
127 final double denominator;
128 if (delta >= 0.0) {
129 // choose a denominator larger in magnitude
130 double dplus = b + JdkMath.sqrt(delta);
131 double dminus = b - JdkMath.sqrt(delta);
132 denominator = JdkMath.abs(dplus) > JdkMath.abs(dminus) ? dplus : dminus;
133 } else {
134 // take the modulus of (B +/- JdkMath.sqrt(delta))
135 denominator = JdkMath.sqrt(b * b - delta);
136 }
137 if (denominator != 0) {
138 x = x2 - 2.0 * c * (x2 - x1) / denominator;
139 // perturb x if it exactly coincides with x1 or x2
140 // the equality tests here are intentional
141 while (x == x1 || x == x2) {
142 x += absoluteAccuracy;
143 }
144 } else {
145 // extremely rare case, get a random number to skip it
146 x = min + JdkMath.random() * (max - min);
147 oldx = Double.POSITIVE_INFINITY;
148 }
149 final double y = computeObjectiveValue(x);
150
151 // check for convergence
152 final double tolerance = JdkMath.max(relativeAccuracy * JdkMath.abs(x), absoluteAccuracy);
153 if (JdkMath.abs(x - oldx) <= tolerance ||
154 JdkMath.abs(y) <= functionValueAccuracy) {
155 return x;
156 }
157
158 // prepare the next iteration
159 x0 = x1;
160 y0 = y1;
161 x1 = x2;
162 y1 = y2;
163 x2 = x;
164 y2 = y;
165 oldx = x;
166 }
167 }
168 }