From fac2e96a318f3ac8b4fe0dc15ac4525165e21af3 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sat, 2 Jul 2011 16:27:03 +0000 Subject: [PATCH] Added BELOW_SIDE and ABOVE_SIDE in the possible allowed solutions for bracketing solvers. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1142244 13f79535-47bb-0310-9956-ffa450edef68 --- .../analysis/solvers/AllowedSolutions.java | 36 ++++++++++++------ .../analysis/solvers/BaseSecantSolver.java | 24 ++++++++---- src/site/xdoc/changes.xml | 3 +- .../solvers/BaseSecantSolverAbstractTest.java | 38 +++++++++++++++++++ 4 files changed, 82 insertions(+), 19 deletions(-) diff --git a/src/main/java/org/apache/commons/math/analysis/solvers/AllowedSolutions.java b/src/main/java/org/apache/commons/math/analysis/solvers/AllowedSolutions.java index 9c729a0ec..ea560907d 100644 --- a/src/main/java/org/apache/commons/math/analysis/solvers/AllowedSolutions.java +++ b/src/main/java/org/apache/commons/math/analysis/solvers/AllowedSolutions.java @@ -18,9 +18,10 @@ package org.apache.commons.math.analysis.solvers; -/** The kinds of solutions that a {@link UnivariateRealSolver (univariate real) - * root-finding algorithm} may accept as solutions. This basically controls - * whether or not under-approximations and over-approximations are allowed. +/** The kinds of solutions that a {@link BracketedUnivariateRealSolver + * (bracketed univariate real) root-finding algorithm} may accept as solutions. + * This basically controls whether or not under-approximations and + * over-approximations are allowed. * *

If all solutions are accepted ({@link #EITHER_SIDE}), then the solution * that the root-finding algorithm returns for a given root may be equal to the @@ -30,8 +31,8 @@ package org.apache.commons.math.analysis.solvers; * tolerances. In certain cases however, in particular for * {@link org.apache.commons.math.ode.events.EventHandler state events} of * {@link org.apache.commons.math.ode.ODEIntegrator ODE solvers}, it - * may be necessary to guarantee that a solution is returned that does not - * under-approximate the solution.

+ * may be necessary to guarantee that a solution is returned that lies on a + * specific side the solution.

* * @see BracketedUnivariateRealSolver * @since 3.0 @@ -40,23 +41,36 @@ package org.apache.commons.math.analysis.solvers; public enum AllowedSolutions { /** There are no additional side restriction on the solutions for * root-finding. That is, both under-approximations and over-approximations - * are allowed. So, if a function f(x) has a root at x = y, then the - * root-finding result s may be smaller than y, equal to y, or greater - * than y. + * are allowed. So, if a function f(x) has a root at x = x0, then the + * root-finding result s may be smaller than x0, equal to x0, or greater + * than x0. */ EITHER_SIDE, /** Only solutions that are less than or equal to the actual root are * acceptable as solutions for root-finding. In other words, * over-approximations are not allowed. So, if a function f(x) has a root - * at x = y, then the root-finding result s must satisfy s <= y. + * at x = x0, then the root-finding result s must satisfy s <= x0. */ LEFT_SIDE, /** Only solutions that are greater than or equal to the actual root are * acceptable as solutions for root-finding. In other words, * under-approximations are not allowed. So, if a function f(x) has a root - * at x = y, then the root-finding result s must satisfy s >= y. + * at x = x0, then the root-finding result s must satisfy s >= x0. */ - RIGHT_SIDE; + RIGHT_SIDE, + + /** Only solutions for which values are less than or equal to zero are + * acceptable as solutions for root-finding. So, if a function f(x) has + * a root at x = x0, then the root-finding result s must satisfy f(s) <= 0. + */ + BELOW_SIDE, + + /** Only solutions for which values are greater than or equal to zero are + * acceptable as solutions for root-finding. So, if a function f(x) has + * a root at x = x0, then the root-finding result s must satisfy f(s) >= 0. + */ + ABOVE_SIDE; + } diff --git a/src/main/java/org/apache/commons/math/analysis/solvers/BaseSecantSolver.java b/src/main/java/org/apache/commons/math/analysis/solvers/BaseSecantSolver.java index cd896b739..3b057919d 100644 --- a/src/main/java/org/apache/commons/math/analysis/solvers/BaseSecantSolver.java +++ b/src/main/java/org/apache/commons/math/analysis/solvers/BaseSecantSolver.java @@ -101,10 +101,6 @@ public abstract class BaseSecantSolver extends AbstractUnivariateRealSolver { final double atol = getAbsoluteAccuracy(); final double rtol = getRelativeAccuracy(); - // Variables to hold new bounds. - double x; - double fx; - // Keep track of inverted intervals, meaning that the left bound is // larger than the right bound. Not used for the original Secant // method. @@ -113,8 +109,8 @@ public abstract class BaseSecantSolver extends AbstractUnivariateRealSolver { // Keep finding better approximations. while (true) { // Calculate the next approximation. - x = x1 - ((f1 * (x1 - x0)) / (f1 - f0)); - fx = computeObjectiveValue(x); + final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0)); + final double fx = computeObjectiveValue(x); // If the new approximation is the exact root, return it. Since // this is not an under-approximation or an over-approximation, @@ -151,7 +147,7 @@ public abstract class BaseSecantSolver extends AbstractUnivariateRealSolver { } // If the function value of the last approximation is too small, - // given the function value accuracy, then we can't get close to + // given the function value accuracy, then we can't get closer to // the root than we already are. if (FastMath.abs(f1) <= ftol) { switch (allowedSolutions) { @@ -167,6 +163,16 @@ public abstract class BaseSecantSolver extends AbstractUnivariateRealSolver { return x1; } break; + case BELOW_SIDE: + if (f1 <= 0) { + return x1; + } + break; + case ABOVE_SIDE: + if (f1 >= 0) { + return x1; + } + break; default: throw new MathInternalError(); } @@ -183,6 +189,10 @@ public abstract class BaseSecantSolver extends AbstractUnivariateRealSolver { return inverted ? x1 : x0; case RIGHT_SIDE: return inverted ? x0 : x1; + case BELOW_SIDE: + return (f1 <= 0) ? x1 : x0; + case ABOVE_SIDE: + return (f1 >= 0) ? x1 : x0; default: throw new MathInternalError(); } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 2b0ba4f2e..bc6030c11 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -59,7 +59,8 @@ The type attribute can be add,update,fix,remove. Modified "SecantSolver" to comply with the original algorithm. Added several - secant-based solvers. + secant-based solvers. Added a way to select the side of the root with bracketing + solvers. Fixed javadoc for ODEIntegrator interface diff --git a/src/test/java/org/apache/commons/math/analysis/solvers/BaseSecantSolverAbstractTest.java b/src/test/java/org/apache/commons/math/analysis/solvers/BaseSecantSolverAbstractTest.java index 3d2199f88..44e3f1793 100644 --- a/src/test/java/org/apache/commons/math/analysis/solvers/BaseSecantSolverAbstractTest.java +++ b/src/test/java/org/apache/commons/math/analysis/solvers/BaseSecantSolverAbstractTest.java @@ -195,4 +195,42 @@ public abstract class BaseSecantSolverAbstractTest { right += 0.3; } } + @Test + public void testSolutionBelowSide() { + UnivariateRealFunction f = new SinFunction(); + UnivariateRealSolver solver = getSolver(); + if (!(solver instanceof BracketedUnivariateRealSolver)) return; + ((BracketedUnivariateRealSolver)solver).setAllowedSolutions(AllowedSolutions.BELOW_SIDE); + double left = -1.5; + double right = 0.05; + for(int i = 0; i < 10; i++) { + // Test whether the allowed solutions are taken into account. + double solution = solver.solve(100, f, left, right); + Assert.assertTrue(f.value(solution) <= 0.0); + + // Prepare for next test. + left -= 0.1; + right += 0.3; + } + } + + @Test + public void testSolutionAboveSide() { + UnivariateRealFunction f = new SinFunction(); + UnivariateRealSolver solver = getSolver(); + if (!(solver instanceof BracketedUnivariateRealSolver)) return; + ((BracketedUnivariateRealSolver)solver).setAllowedSolutions(AllowedSolutions.ABOVE_SIDE); + double left = -1.5; + double right = 0.05; + for(int i = 0; i < 10; i++) { + // Test whether the allowed solutions are taken into account. + double solution = solver.solve(100, f, left, right); + Assert.assertTrue(f.value(solution) >= 0.0); + + // Prepare for next test. + left -= 0.1; + right += 0.3; + } + } + }