Repository: commons-math Updated Branches: refs/heads/feature-MATH-1333 [created] 160696e7f
MATH-1333 Unit test showing the problem. Thanks to Connor Petty for the report. Assumption made in the code is wrong. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/160696e7 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/160696e7 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/160696e7 Branch: refs/heads/feature-MATH-1333 Commit: 160696e7faa591e14658aec12fc159efdb65cf12 Parents: 12c9a04 Author: Gilles <er...@apache.org> Authored: Wed Mar 9 16:47:09 2016 +0100 Committer: Gilles <er...@apache.org> Committed: Wed Mar 9 16:47:09 2016 +0100 ---------------------------------------------------------------------- .../math4/analysis/solvers/MullerSolver.java | 9 ++++-- .../analysis/solvers/MullerSolverTest.java | 31 ++++++++++++++++++++ 2 files changed, 38 insertions(+), 2 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/160696e7/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java b/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java index b714399..d7fc189 100644 --- a/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java +++ b/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java @@ -29,8 +29,6 @@ import org.apache.commons.math4.util.FastMath; * <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, @@ -161,6 +159,13 @@ public class MullerSolver extends AbstractUnivariateSolver { // 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; + + // XXX debug + if (!isSequence(x0, x, x2)) { + System.out.println("x=" + x + " x0=" + x0 + " x2=" + x2); + throw new org.apache.commons.math4.exception.MathInternalError(); + } + final double y = computeObjectiveValue(x); // check for convergence http://git-wip-us.apache.org/repos/asf/commons-math/blob/160696e7/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java b/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java index e82be52..2bb0b1b 100644 --- a/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java +++ b/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java @@ -147,4 +147,35 @@ public final class MullerSolverTest { // expected } } + + @Test + public void testMath1333() { + final UnivariateFunction logFunction = new UnivariateFunction() { + private double log1pe(double x) { + if (x > 0) { + return x + FastMath.log1p(FastMath.exp(-x)); + } else { + return FastMath.log1p(FastMath.exp(x)); + } + } + + @Override + public double value(double x) { + final double a = 0.15076136473214652; + final double b = 4.880819340168248; + final double c = -2330.4196672490493; + final double d = 1.1871451743330544E-16; + //aa*log(1+e^(bbx+c))+d - 0.01 * x - 20 * 0.01 + return a * a * log1pe(b * b * x + c) + d - 0.01 * x - 20 * 0.01; + } + }; + + final UnivariateSolver solver = new MullerSolver(0.25); + final double min = 20; + final double max = 100.04173804515072; + final double result = solver.solve(1000, logFunction, min, max, 100 / (double) 3); + + Assert.assertTrue(result + " < " + min, result >= min); + Assert.assertTrue(result + " > " + max, result <= max); + } }
