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