diff options
author | Karl Shaffer <karlshaffer@google.com> | 2023-08-10 20:17:10 +0000 |
---|---|---|
committer | Automerger Merge Worker <android-build-automerger-merge-worker@system.gserviceaccount.com> | 2023-08-10 20:17:10 +0000 |
commit | 0018f64b87f6b2df2a4c866aa3eb8ab2e45efbea (patch) | |
tree | ace24ba4307d4978ee3134f7da671a77ad172da0 /src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java | |
parent | 21944074de3a1f7bafb37541992e7277d639fda3 (diff) | |
parent | 1354beaf452fc42c23c82fc80d2f2e9ffcbf3688 (diff) | |
download | apache-commons-math-0018f64b87f6b2df2a4c866aa3eb8ab2e45efbea.tar.gz |
Check-in commons-math 3.6.1 am: 1354beaf45
Original change: https://android-review.googlesource.com/c/platform/external/apache-commons-math/+/2702413
Change-Id: I4990b1b31e5dd9596f208d93dd22cc8d5e470689
Signed-off-by: Automerger Merge Worker <android-build-automerger-merge-worker@system.gserviceaccount.com>
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java')
-rw-r--r-- | src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java | 202 |
1 files changed, 202 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java new file mode 100644 index 0000000..2e59ed4 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver.java @@ -0,0 +1,202 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.solvers; + +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; + +/** + * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html"> + * Muller's Method</a> for root finding of real univariate functions. For + * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477, + * chapter 3. + * <p> + * Muller's method applies to both real and complex functions, but here we + * restrict ourselves to real functions. + * This class differs from {@link MullerSolver} in the way it avoids complex + * operations.</p><p> + * Muller's original method would have function evaluation at complex point. + * Since our f(x) is real, we have to find ways to avoid that. Bracketing + * condition is one way to go: by requiring bracketing in every iteration, + * the newly computed approximation is guaranteed to be real.</p> + * <p> + * Normally Muller's method converges quadratically in the vicinity of a + * zero, however it may be very slow in regions far away from zeros. For + * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use + * bisection as a safety backup if it performs very poorly.</p> + * <p> + * The formulas here use divided differences directly.</p> + * + * @since 1.2 + * @see MullerSolver2 + */ +public class MullerSolver extends AbstractUnivariateSolver { + + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver with default accuracy (1e-6). + */ + public MullerSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public MullerSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public MullerSolver(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException, + NumberIsTooLargeException, + NoBracketingException { + final double min = getMin(); + final double max = getMax(); + final double initial = getStartValue(); + + final double functionValueAccuracy = getFunctionValueAccuracy(); + + verifySequence(min, initial, max); + + // check for zeros before verifying bracketing + final double fMin = computeObjectiveValue(min); + if (FastMath.abs(fMin) < functionValueAccuracy) { + return min; + } + final double fMax = computeObjectiveValue(max); + if (FastMath.abs(fMax) < functionValueAccuracy) { + return max; + } + final double fInitial = computeObjectiveValue(initial); + if (FastMath.abs(fInitial) < functionValueAccuracy) { + return initial; + } + + verifyBracketing(min, max); + + if (isBracketing(min, initial)) { + return solve(min, initial, fMin, fInitial); + } else { + return solve(initial, max, fInitial, fMax); + } + } + + /** + * Find a real root in the given interval. + * + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param fMin function value at the lower bound. + * @param fMax function value at the upper bound. + * @return the point at which the function value is zero. + * @throws TooManyEvaluationsException if the allowed number of calls to + * the function to be solved has been exhausted. + */ + private double solve(double min, double max, + double fMin, double fMax) + throws TooManyEvaluationsException { + final double relativeAccuracy = getRelativeAccuracy(); + final double absoluteAccuracy = getAbsoluteAccuracy(); + final double functionValueAccuracy = getFunctionValueAccuracy(); + + // [x0, x2] is the bracketing interval in each iteration + // x1 is the last approximation and an interpolation point in (x0, x2) + // x is the new root approximation and new x1 for next round + // d01, d12, d012 are divided differences + + double x0 = min; + double y0 = fMin; + double x2 = max; + double y2 = fMax; + double x1 = 0.5 * (x0 + x2); + double y1 = computeObjectiveValue(x1); + + double oldx = Double.POSITIVE_INFINITY; + while (true) { + // Muller's method employs quadratic interpolation through + // x0, x1, x2 and x is the zero of the interpolating parabola. + // Due to bracketing condition, this parabola must have two + // real roots and we choose one in [x0, x2] to be x. + final double d01 = (y1 - y0) / (x1 - x0); + final double d12 = (y2 - y1) / (x2 - x1); + final double d012 = (d12 - d01) / (x2 - x0); + final double c1 = d01 + (x1 - x0) * d012; + final double delta = c1 * c1 - 4 * y1 * d012; + final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta)); + final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta)); + // xplus and xminus are two roots of parabola and at least + // one of them should lie in (x0, x2) + final double x = isSequence(x0, xplus, x2) ? xplus : xminus; + final double y = computeObjectiveValue(x); + + // check for convergence + final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); + if (FastMath.abs(x - oldx) <= tolerance || + FastMath.abs(y) <= functionValueAccuracy) { + return x; + } + + // Bisect if convergence is too slow. Bisection would waste + // our calculation of x, hopefully it won't happen often. + // the real number equality test x == x1 is intentional and + // completes the proximity tests above it + boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) || + (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) || + (x == x1); + // prepare the new bracketing interval for next iteration + if (!bisect) { + x0 = x < x1 ? x0 : x1; + y0 = x < x1 ? y0 : y1; + x2 = x > x1 ? x2 : x1; + y2 = x > x1 ? y2 : y1; + x1 = x; y1 = y; + oldx = x; + } else { + double xm = 0.5 * (x0 + x2); + double ym = computeObjectiveValue(xm); + if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) { + x2 = xm; y2 = ym; + } else { + x0 = xm; y0 = ym; + } + x1 = 0.5 * (x0 + x2); + y1 = computeObjectiveValue(x1); + oldx = Double.POSITIVE_INFINITY; + } + } + } +} |