diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/solvers')
30 files changed, 4677 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractDifferentiableUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractDifferentiableUnivariateSolver.java new file mode 100644 index 0000000..d0fda00 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractDifferentiableUnivariateSolver.java @@ -0,0 +1,82 @@ +/* + * 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.analysis.DifferentiableUnivariateFunction; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Provide a default implementation for several functions useful to generic + * solvers. + * + * @since 3.0 + * @deprecated as of 3.1, replaced by {@link AbstractUnivariateDifferentiableSolver} + */ +@Deprecated +public abstract class AbstractDifferentiableUnivariateSolver + extends BaseAbstractUnivariateSolver<DifferentiableUnivariateFunction> + implements DifferentiableUnivariateSolver { + /** Derivative of the function to solve. */ + private UnivariateFunction functionDerivative; + + /** + * Construct a solver with given absolute accuracy. + * + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractDifferentiableUnivariateSolver(final double absoluteAccuracy) { + super(absoluteAccuracy); + } + + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + */ + protected AbstractDifferentiableUnivariateSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } + + /** + * Compute the objective function value. + * + * @param point Point at which the objective function must be evaluated. + * @return the objective function value at specified point. + * @throws TooManyEvaluationsException if the maximal number of evaluations is exceeded. + */ + protected double computeDerivativeObjectiveValue(double point) + throws TooManyEvaluationsException { + incrementEvaluationCount(); + return functionDerivative.value(point); + } + + /** + * {@inheritDoc} + */ + @Override + protected void setup(int maxEval, DifferentiableUnivariateFunction f, + double min, double max, double startValue) { + super.setup(maxEval, f, min, max, startValue); + functionDerivative = f.derivative(); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractPolynomialSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractPolynomialSolver.java new file mode 100644 index 0000000..d641e87 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractPolynomialSolver.java @@ -0,0 +1,80 @@ +/* + * 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.analysis.polynomials.PolynomialFunction; + +/** + * Base class for solvers. + * + * @since 3.0 + */ +public abstract class AbstractPolynomialSolver + extends BaseAbstractUnivariateSolver<PolynomialFunction> + implements PolynomialSolver { + /** Function. */ + private PolynomialFunction polynomialFunction; + + /** + * Construct a solver with given absolute accuracy. + * + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractPolynomialSolver(final double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractPolynomialSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + */ + protected AbstractPolynomialSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected void setup(int maxEval, PolynomialFunction f, + double min, double max, double startValue) { + super.setup(maxEval, f, min, max, startValue); + polynomialFunction = f; + } + + /** + * @return the coefficients of the polynomial function. + */ + protected double[] getCoefficients() { + return polynomialFunction.getCoefficients(); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateDifferentiableSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateDifferentiableSolver.java new file mode 100644 index 0000000..9745e9b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateDifferentiableSolver.java @@ -0,0 +1,82 @@ +/* + * 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.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Provide a default implementation for several functions useful to generic + * solvers. + * + * @since 3.1 + */ +public abstract class AbstractUnivariateDifferentiableSolver + extends BaseAbstractUnivariateSolver<UnivariateDifferentiableFunction> + implements UnivariateDifferentiableSolver { + + /** Function to solve. */ + private UnivariateDifferentiableFunction function; + + /** + * Construct a solver with given absolute accuracy. + * + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractUnivariateDifferentiableSolver(final double absoluteAccuracy) { + super(absoluteAccuracy); + } + + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + */ + protected AbstractUnivariateDifferentiableSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } + + /** + * Compute the objective function value. + * + * @param point Point at which the objective function must be evaluated. + * @return the objective function value and derivative at specified point. + * @throws TooManyEvaluationsException + * if the maximal number of evaluations is exceeded. + */ + protected DerivativeStructure computeObjectiveValueAndDerivative(double point) + throws TooManyEvaluationsException { + incrementEvaluationCount(); + return function.value(new DerivativeStructure(1, 1, 0, point)); + } + + /** + * {@inheritDoc} + */ + @Override + protected void setup(int maxEval, UnivariateDifferentiableFunction f, + double min, double max, double startValue) { + super.setup(maxEval, f, min, max, startValue); + function = f; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateSolver.java new file mode 100644 index 0000000..078c70f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/AbstractUnivariateSolver.java @@ -0,0 +1,60 @@ +/* + * 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.analysis.UnivariateFunction; + +/** + * Base class for solvers. + * + * @since 3.0 + */ +public abstract class AbstractUnivariateSolver + extends BaseAbstractUnivariateSolver<UnivariateFunction> + implements UnivariateSolver { + /** + * Construct a solver with given absolute accuracy. + * + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractUnivariateSolver(final double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + */ + protected AbstractUnivariateSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + */ + protected AbstractUnivariateSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/AllowedSolution.java b/src/main/java/org/apache/commons/math3/analysis/solvers/AllowedSolution.java new file mode 100644 index 0000000..a02a29b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/AllowedSolution.java @@ -0,0 +1,75 @@ +/* + * 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; + + +/** The kinds of solutions that a {@link BracketedUnivariateSolver + * (bracketed univariate real) root-finding algorithm} may accept as solutions. + * This basically controls whether or not under-approximations and + * over-approximations are allowed. + * + * <p>If all solutions are accepted ({@link #ANY_SIDE}), then the solution + * that the root-finding algorithm returns for a given root may be equal to the + * actual root, but it may also be an approximation that is slightly smaller + * or slightly larger than the actual root. Root-finding algorithms generally + * only guarantee that the returned solution is within the requested + * tolerances. In certain cases however, in particular for + * {@link org.apache.commons.math3.ode.events.EventHandler state events} of + * {@link org.apache.commons.math3.ode.ODEIntegrator ODE solvers}, it + * may be necessary to guarantee that a solution is returned that lies on a + * specific side the solution.</p> + * + * @see BracketedUnivariateSolver + * @since 3.0 + */ +public enum AllowedSolution { + /** 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 = x0, then the + * root-finding result s may be smaller than x0, equal to x0, or greater + * than x0. + */ + ANY_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 = 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 = x0, then the root-finding result s must satisfy s >= x0. + */ + 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/math3/analysis/solvers/BaseAbstractUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseAbstractUnivariateSolver.java new file mode 100644 index 0000000..12b30c6 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseAbstractUnivariateSolver.java @@ -0,0 +1,318 @@ +/* + * 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.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.util.IntegerSequence; +import org.apache.commons.math3.util.MathUtils; + +/** + * Provide a default implementation for several functions useful to generic + * solvers. + * The default values for relative and function tolerances are 1e-14 + * and 1e-15, respectively. It is however highly recommended to not + * rely on the default, but rather carefully consider values that match + * user's expectations, as well as the specifics of each implementation. + * + * @param <FUNC> Type of function to solve. + * + * @since 2.0 + */ +public abstract class BaseAbstractUnivariateSolver<FUNC extends UnivariateFunction> + implements BaseUnivariateSolver<FUNC> { + /** Default relative accuracy. */ + private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14; + /** Default function value accuracy. */ + private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15; + /** Function value accuracy. */ + private final double functionValueAccuracy; + /** Absolute accuracy. */ + private final double absoluteAccuracy; + /** Relative accuracy. */ + private final double relativeAccuracy; + /** Evaluations counter. */ + private IntegerSequence.Incrementor evaluations; + /** Lower end of search interval. */ + private double searchMin; + /** Higher end of search interval. */ + private double searchMax; + /** Initial guess. */ + private double searchStart; + /** Function to solve. */ + private FUNC function; + + /** + * Construct a solver with given absolute accuracy. + * + * @param absoluteAccuracy Maximum absolute error. + */ + protected BaseAbstractUnivariateSolver(final double absoluteAccuracy) { + this(DEFAULT_RELATIVE_ACCURACY, + absoluteAccuracy, + DEFAULT_FUNCTION_VALUE_ACCURACY); + } + + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + */ + protected BaseAbstractUnivariateSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + this(relativeAccuracy, + absoluteAccuracy, + DEFAULT_FUNCTION_VALUE_ACCURACY); + } + + /** + * Construct a solver with given accuracies. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + */ + protected BaseAbstractUnivariateSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + this.absoluteAccuracy = absoluteAccuracy; + this.relativeAccuracy = relativeAccuracy; + this.functionValueAccuracy = functionValueAccuracy; + this.evaluations = IntegerSequence.Incrementor.create(); + } + + /** {@inheritDoc} */ + public int getMaxEvaluations() { + return evaluations.getMaximalCount(); + } + /** {@inheritDoc} */ + public int getEvaluations() { + return evaluations.getCount(); + } + /** + * @return the lower end of the search interval. + */ + public double getMin() { + return searchMin; + } + /** + * @return the higher end of the search interval. + */ + public double getMax() { + return searchMax; + } + /** + * @return the initial guess. + */ + public double getStartValue() { + return searchStart; + } + /** + * {@inheritDoc} + */ + public double getAbsoluteAccuracy() { + return absoluteAccuracy; + } + /** + * {@inheritDoc} + */ + public double getRelativeAccuracy() { + return relativeAccuracy; + } + /** + * {@inheritDoc} + */ + public double getFunctionValueAccuracy() { + return functionValueAccuracy; + } + + /** + * Compute the objective function value. + * + * @param point Point at which the objective function must be evaluated. + * @return the objective function value at specified point. + * @throws TooManyEvaluationsException if the maximal number of evaluations + * is exceeded. + */ + protected double computeObjectiveValue(double point) + throws TooManyEvaluationsException { + incrementEvaluationCount(); + return function.value(point); + } + + /** + * Prepare for computation. + * Subclasses must call this method if they override any of the + * {@code solve} methods. + * + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @param maxEval Maximum number of evaluations. + * @exception NullArgumentException if f is null + */ + protected void setup(int maxEval, + FUNC f, + double min, double max, + double startValue) + throws NullArgumentException { + // Checks. + MathUtils.checkNotNull(f); + + // Reset. + searchMin = min; + searchMax = max; + searchStart = startValue; + function = f; + evaluations = evaluations.withMaximalCount(maxEval).withStart(0); + } + + /** {@inheritDoc} */ + public double solve(int maxEval, FUNC f, double min, double max, double startValue) + throws TooManyEvaluationsException, + NoBracketingException { + // Initialization. + setup(maxEval, f, min, max, startValue); + + // Perform computation. + return doSolve(); + } + + /** {@inheritDoc} */ + public double solve(int maxEval, FUNC f, double min, double max) { + return solve(maxEval, f, min, max, min + 0.5 * (max - min)); + } + + /** {@inheritDoc} */ + public double solve(int maxEval, FUNC f, double startValue) + throws TooManyEvaluationsException, + NoBracketingException { + return solve(maxEval, f, Double.NaN, Double.NaN, startValue); + } + + /** + * Method for implementing actual optimization algorithms in derived + * classes. + * + * @return the root. + * @throws TooManyEvaluationsException if the maximal number of evaluations + * is exceeded. + * @throws NoBracketingException if the initial search interval does not bracket + * a root and the solver requires it. + */ + protected abstract double doSolve() + throws TooManyEvaluationsException, NoBracketingException; + + /** + * Check whether the function takes opposite signs at the endpoints. + * + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @return {@code true} if the function values have opposite signs at the + * given points. + */ + protected boolean isBracketing(final double lower, + final double upper) { + return UnivariateSolverUtils.isBracketing(function, lower, upper); + } + + /** + * Check whether the arguments form a (strictly) increasing sequence. + * + * @param start First number. + * @param mid Second number. + * @param end Third number. + * @return {@code true} if the arguments form an increasing sequence. + */ + protected boolean isSequence(final double start, + final double mid, + final double end) { + return UnivariateSolverUtils.isSequence(start, mid, end); + } + + /** + * Check that the endpoints specify an interval. + * + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + protected void verifyInterval(final double lower, + final double upper) + throws NumberIsTooLargeException { + UnivariateSolverUtils.verifyInterval(lower, upper); + } + + /** + * Check that {@code lower < initial < upper}. + * + * @param lower Lower endpoint. + * @param initial Initial value. + * @param upper Upper endpoint. + * @throws NumberIsTooLargeException if {@code lower >= initial} or + * {@code initial >= upper}. + */ + protected void verifySequence(final double lower, + final double initial, + final double upper) + throws NumberIsTooLargeException { + UnivariateSolverUtils.verifySequence(lower, initial, upper); + } + + /** + * Check that the endpoints specify an interval and the function takes + * opposite signs at the endpoints. + * + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @throws NullArgumentException if the function has not been set. + * @throws NoBracketingException if the function has the same sign at + * the endpoints. + */ + protected void verifyBracketing(final double lower, + final double upper) + throws NullArgumentException, + NoBracketingException { + UnivariateSolverUtils.verifyBracketing(function, lower, upper); + } + + /** + * Increment the evaluation count by one. + * Method {@link #computeObjectiveValue(double)} calls this method internally. + * It is provided for subclasses that do not exclusively use + * {@code computeObjectiveValue} to solve the function. + * See e.g. {@link AbstractUnivariateDifferentiableSolver}. + * + * @throws TooManyEvaluationsException when the allowed number of function + * evaluations has been exhausted. + */ + protected void incrementEvaluationCount() + throws TooManyEvaluationsException { + try { + evaluations.increment(); + } catch (MaxCountExceededException e) { + throw new TooManyEvaluationsException(e.getMax()); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BaseSecantSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseSecantSolver.java new file mode 100644 index 0000000..44a2173 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseSecantSolver.java @@ -0,0 +1,278 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.MathInternalError; + +/** + * Base class for all bracketing <em>Secant</em>-based methods for root-finding + * (approximating a zero of a univariate real function). + * + * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and + * {@link IllinoisSolver <em>Illinois</em>} methods is based on the + * following article: M. Dowell and P. Jarratt, + * <em>A modified regula falsi method for computing the root of an + * equation</em>, BIT Numerical Mathematics, volume 11, number 2, + * pages 168-174, Springer, 1971.</p> + * + * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is + * based on the following article: M. Dowell and P. Jarratt, + * <em>The "Pegasus" method for computing the root of an equation</em>, + * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer, + * 1972.</p> + * + * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a + * bracketing method, so it is not implemented here. It has a separate + * implementation.</p> + * + * @since 3.0 + */ +public abstract class BaseSecantSolver + extends AbstractUnivariateSolver + implements BracketedUnivariateSolver<UnivariateFunction> { + + /** Default absolute accuracy. */ + protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** The kinds of solutions that the algorithm may accept. */ + private AllowedSolution allowed; + + /** The <em>Secant</em>-based root-finding method to use. */ + private final Method method; + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + * @param method <em>Secant</em>-based root-finding method to use. + */ + protected BaseSecantSolver(final double absoluteAccuracy, final Method method) { + super(absoluteAccuracy); + this.allowed = AllowedSolution.ANY_SIDE; + this.method = method; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param method <em>Secant</em>-based root-finding method to use. + */ + protected BaseSecantSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final Method method) { + super(relativeAccuracy, absoluteAccuracy); + this.allowed = AllowedSolution.ANY_SIDE; + this.method = method; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Maximum relative error. + * @param absoluteAccuracy Maximum absolute error. + * @param functionValueAccuracy Maximum function value error. + * @param method <em>Secant</em>-based root-finding method to use + */ + protected BaseSecantSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy, + final Method method) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + this.allowed = AllowedSolution.ANY_SIDE; + this.method = method; + } + + /** {@inheritDoc} */ + public double solve(final int maxEval, final UnivariateFunction f, + final double min, final double max, + final AllowedSolution allowedSolution) { + return solve(maxEval, f, min, max, min + 0.5 * (max - min), allowedSolution); + } + + /** {@inheritDoc} */ + public double solve(final int maxEval, final UnivariateFunction f, + final double min, final double max, final double startValue, + final AllowedSolution allowedSolution) { + this.allowed = allowedSolution; + return super.solve(maxEval, f, min, max, startValue); + } + + /** {@inheritDoc} */ + @Override + public double solve(final int maxEval, final UnivariateFunction f, + final double min, final double max, final double startValue) { + return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE); + } + + /** + * {@inheritDoc} + * + * @throws ConvergenceException if the algorithm failed due to finite + * precision. + */ + @Override + protected final double doSolve() + throws ConvergenceException { + // Get initial solution + double x0 = getMin(); + double x1 = getMax(); + double f0 = computeObjectiveValue(x0); + double f1 = computeObjectiveValue(x1); + + // If one of the bounds is the exact root, return it. Since these are + // not under-approximations or over-approximations, we can return them + // regardless of the allowed solutions. + if (f0 == 0.0) { + return x0; + } + if (f1 == 0.0) { + return x1; + } + + // Verify bracketing of initial solution. + verifyBracketing(x0, x1); + + // Get accuracies. + final double ftol = getFunctionValueAccuracy(); + final double atol = getAbsoluteAccuracy(); + final double rtol = getRelativeAccuracy(); + + // Keep track of inverted intervals, meaning that the left bound is + // larger than the right bound. + boolean inverted = false; + + // Keep finding better approximations. + while (true) { + // Calculate the next approximation. + 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, + // we can return it regardless of the allowed solutions. + if (fx == 0.0) { + return x; + } + + // Update the bounds with the new approximation. + if (f1 * fx < 0) { + // The value of x1 has switched to the other bound, thus inverting + // the interval. + x0 = x1; + f0 = f1; + inverted = !inverted; + } else { + switch (method) { + case ILLINOIS: + f0 *= 0.5; + break; + case PEGASUS: + f0 *= f1 / (f1 + fx); + break; + case REGULA_FALSI: + // Detect early that algorithm is stuck, instead of waiting + // for the maximum number of iterations to be exceeded. + if (x == x1) { + throw new ConvergenceException(); + } + break; + default: + // Should never happen. + throw new MathInternalError(); + } + } + // Update from [x0, x1] to [x0, x]. + x1 = x; + f1 = fx; + + // If the function value of the last approximation is too small, + // given the function value accuracy, then we can't get closer to + // the root than we already are. + if (FastMath.abs(f1) <= ftol) { + switch (allowed) { + case ANY_SIDE: + return x1; + case LEFT_SIDE: + if (inverted) { + return x1; + } + break; + case RIGHT_SIDE: + if (!inverted) { + 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(); + } + } + + // If the current interval is within the given accuracies, we + // are satisfied with the current approximation. + if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), + atol)) { + switch (allowed) { + case ANY_SIDE: + return x1; + case LEFT_SIDE: + 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(); + } + } + } + } + + /** <em>Secant</em>-based root-finding methods. */ + protected enum Method { + + /** + * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or + * <em>False Position</em> method. + */ + REGULA_FALSI, + + /** The {@link IllinoisSolver <em>Illinois</em>} method. */ + ILLINOIS, + + /** The {@link PegasusSolver <em>Pegasus</em>} method. */ + PEGASUS; + + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BaseUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseUnivariateSolver.java new file mode 100644 index 0000000..f00590e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BaseUnivariateSolver.java @@ -0,0 +1,142 @@ +/* + * 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.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + + +/** + * Interface for (univariate real) rootfinding algorithms. + * Implementations will search for only one zero in the given interval. + * + * This class is not intended for use outside of the Apache Commons Math + * library, regular user should rely on more specific interfaces like + * {@link UnivariateSolver}, {@link PolynomialSolver} or {@link + * DifferentiableUnivariateSolver}. + * @param <FUNC> Type of function to solve. + * + * @since 3.0 + * @see UnivariateSolver + * @see PolynomialSolver + * @see DifferentiableUnivariateSolver + */ +public interface BaseUnivariateSolver<FUNC extends UnivariateFunction> { + /** + * Get the maximum number of function evaluations. + * + * @return the maximum number of function evaluations. + */ + int getMaxEvaluations(); + + /** + * Get the number of evaluations of the objective function. + * The number of evaluations corresponds to the last call to the + * {@code optimize} method. It is 0 if the method has not been + * called yet. + * + * @return the number of evaluations of the objective function. + */ + int getEvaluations(); + + /** + * Get the absolute accuracy of the solver. Solutions returned by the + * solver should be accurate to this tolerance, i.e., if ε is the + * absolute accuracy of the solver and {@code v} is a value returned by + * one of the {@code solve} methods, then a root of the function should + * exist somewhere in the interval ({@code v} - ε, {@code v} + ε). + * + * @return the absolute accuracy. + */ + double getAbsoluteAccuracy(); + + /** + * Get the relative accuracy of the solver. The contract for relative + * accuracy is the same as {@link #getAbsoluteAccuracy()}, but using + * relative, rather than absolute error. If ρ is the relative accuracy + * configured for a solver and {@code v} is a value returned, then a root + * of the function should exist somewhere in the interval + * ({@code v} - ρ {@code v}, {@code v} + ρ {@code v}). + * + * @return the relative accuracy. + */ + double getRelativeAccuracy(); + + /** + * Get the function value accuracy of the solver. If {@code v} is + * a value returned by the solver for a function {@code f}, + * then by contract, {@code |f(v)|} should be less than or equal to + * the function value accuracy configured for the solver. + * + * @return the function value accuracy. + */ + double getFunctionValueAccuracy(); + + /** + * Solve for a zero root in the given interval. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @return a value where the function is zero. + * @throws MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + double solve(int maxEval, FUNC f, double min, double max) + throws MathIllegalArgumentException, TooManyEvaluationsException; + + /** + * Solve for a zero in the given interval, start at {@code startValue}. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @return a value where the function is zero. + * @throws MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + double solve(int maxEval, FUNC f, double min, double max, double startValue) + throws MathIllegalArgumentException, TooManyEvaluationsException; + + /** + * Solve for a zero in the vicinity of {@code startValue}. + * + * @param f Function to solve. + * @param startValue Start value to use. + * @return a value where the function is zero. + * @param maxEval Maximum number of evaluations. + * @throws org.apache.commons.math3.exception.MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + double solve(int maxEval, FUNC f, double startValue); +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BisectionSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BisectionSolver.java new file mode 100644 index 0000000..49f4057 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BisectionSolver.java @@ -0,0 +1,91 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Implements the <a href="http://mathworld.wolfram.com/Bisection.html"> + * bisection algorithm</a> for finding zeros of univariate real functions. + * <p> + * The function should be continuous but not necessarily smooth.</p> + * + */ +public class BisectionSolver extends AbstractUnivariateSolver { + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver with default accuracy (1e-6). + */ + public BisectionSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public BisectionSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public BisectionSolver(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException { + double min = getMin(); + double max = getMax(); + verifyInterval(min, max); + final double absoluteAccuracy = getAbsoluteAccuracy(); + double m; + double fm; + double fmin; + + while (true) { + m = UnivariateSolverUtils.midpoint(min, max); + fmin = computeObjectiveValue(min); + fm = computeObjectiveValue(m); + + if (fm * fmin > 0) { + // max and m bracket the root. + min = m; + } else { + // min and m bracket the root. + max = m; + } + + if (FastMath.abs(max - min) <= absoluteAccuracy) { + m = UnivariateSolverUtils.midpoint(min, max); + return m; + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedRealFieldUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedRealFieldUnivariateSolver.java new file mode 100644 index 0000000..55562e5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedRealFieldUnivariateSolver.java @@ -0,0 +1,142 @@ +/* + * 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.RealFieldElement; +import org.apache.commons.math3.analysis.RealFieldUnivariateFunction; + +/** Interface for {@link UnivariateSolver (univariate real) root-finding + * algorithms} that maintain a bracketed solution. There are several advantages + * to having such root-finding algorithms: + * <ul> + * <li>The bracketed solution guarantees that the root is kept within the + * interval. As such, these algorithms generally also guarantee + * convergence.</li> + * <li>The bracketed solution means that we have the opportunity to only + * return roots that are greater than or equal to the actual root, or + * are less than or equal to the actual root. That is, we can control + * whether under-approximations and over-approximations are + * {@link AllowedSolution allowed solutions}. Other root-finding + * algorithms can usually only guarantee that the solution (the root that + * was found) is around the actual root.</li> + * </ul> + * + * <p>For backwards compatibility, all root-finding algorithms must have + * {@link AllowedSolution#ANY_SIDE ANY_SIDE} as default for the allowed + * solutions.</p> + * + * @see AllowedSolution + * @param <T> the type of the field elements + * @since 3.6 + */ +public interface BracketedRealFieldUnivariateSolver<T extends RealFieldElement<T>> { + + /** + * Get the maximum number of function evaluations. + * + * @return the maximum number of function evaluations. + */ + int getMaxEvaluations(); + + /** + * Get the number of evaluations of the objective function. + * The number of evaluations corresponds to the last call to the + * {@code optimize} method. It is 0 if the method has not been + * called yet. + * + * @return the number of evaluations of the objective function. + */ + int getEvaluations(); + + /** + * Get the absolute accuracy of the solver. Solutions returned by the + * solver should be accurate to this tolerance, i.e., if ε is the + * absolute accuracy of the solver and {@code v} is a value returned by + * one of the {@code solve} methods, then a root of the function should + * exist somewhere in the interval ({@code v} - ε, {@code v} + ε). + * + * @return the absolute accuracy. + */ + T getAbsoluteAccuracy(); + + /** + * Get the relative accuracy of the solver. The contract for relative + * accuracy is the same as {@link #getAbsoluteAccuracy()}, but using + * relative, rather than absolute error. If ρ is the relative accuracy + * configured for a solver and {@code v} is a value returned, then a root + * of the function should exist somewhere in the interval + * ({@code v} - ρ {@code v}, {@code v} + ρ {@code v}). + * + * @return the relative accuracy. + */ + T getRelativeAccuracy(); + + /** + * Get the function value accuracy of the solver. If {@code v} is + * a value returned by the solver for a function {@code f}, + * then by contract, {@code |f(v)|} should be less than or equal to + * the function value accuracy configured for the solver. + * + * @return the function value accuracy. + */ + T getFunctionValueAccuracy(); + + /** + * Solve for a zero in the given interval. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return A value where the function is zero. + * @throws org.apache.commons.math3.exception.MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + T solve(int maxEval, RealFieldUnivariateFunction<T> f, T min, T max, + AllowedSolution allowedSolution); + + /** + * Solve for a zero in the given interval, start at {@code startValue}. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return A value where the function is zero. + * @throws org.apache.commons.math3.exception.MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + T solve(int maxEval, RealFieldUnivariateFunction<T> f, T min, T max, T startValue, + AllowedSolution allowedSolution); + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedUnivariateSolver.java new file mode 100644 index 0000000..789fc99 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketedUnivariateSolver.java @@ -0,0 +1,92 @@ +/* + * 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.analysis.UnivariateFunction; + +/** Interface for {@link UnivariateSolver (univariate real) root-finding + * algorithms} that maintain a bracketed solution. There are several advantages + * to having such root-finding algorithms: + * <ul> + * <li>The bracketed solution guarantees that the root is kept within the + * interval. As such, these algorithms generally also guarantee + * convergence.</li> + * <li>The bracketed solution means that we have the opportunity to only + * return roots that are greater than or equal to the actual root, or + * are less than or equal to the actual root. That is, we can control + * whether under-approximations and over-approximations are + * {@link AllowedSolution allowed solutions}. Other root-finding + * algorithms can usually only guarantee that the solution (the root that + * was found) is around the actual root.</li> + * </ul> + * + * <p>For backwards compatibility, all root-finding algorithms must have + * {@link AllowedSolution#ANY_SIDE ANY_SIDE} as default for the allowed + * solutions.</p> + * @param <FUNC> Type of function to solve. + * + * @see AllowedSolution + * @since 3.0 + */ +public interface BracketedUnivariateSolver<FUNC extends UnivariateFunction> + extends BaseUnivariateSolver<FUNC> { + + /** + * Solve for a zero in the given interval. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return A value where the function is zero. + * @throws org.apache.commons.math3.exception.MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + double solve(int maxEval, FUNC f, double min, double max, + AllowedSolution allowedSolution); + + /** + * Solve for a zero in the given interval, start at {@code startValue}. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return A value where the function is zero. + * @throws org.apache.commons.math3.exception.MathIllegalArgumentException + * if the arguments do not satisfy the requirements specified by the solver. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if + * the allowed number of evaluations is exceeded. + */ + double solve(int maxEval, FUNC f, double min, double max, double startValue, + AllowedSolution allowedSolution); + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java new file mode 100644 index 0000000..956bf9e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver.java @@ -0,0 +1,411 @@ +/* + * 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.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; + +/** + * This class implements a modification of the <a + * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>. + * <p> + * The changes with respect to the original Brent algorithm are: + * <ul> + * <li>the returned value is chosen in the current interval according + * to user specified {@link AllowedSolution},</li> + * <li>the maximal order for the invert polynomial root search is + * user-specified instead of being invert quadratic only</li> + * </ul><p> + * The given interval must bracket the root.</p> + * + */ +public class BracketingNthOrderBrentSolver + extends AbstractUnivariateSolver + implements BracketedUnivariateSolver<UnivariateFunction> { + + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** Default maximal order. */ + private static final int DEFAULT_MAXIMAL_ORDER = 5; + + /** Maximal aging triggering an attempt to balance the bracketing interval. */ + private static final int MAXIMAL_AGING = 2; + + /** Reduction factor for attempts to balance the bracketing interval. */ + private static final double REDUCTION_FACTOR = 1.0 / 16.0; + + /** Maximal order. */ + private final int maximalOrder; + + /** The kinds of solutions that the algorithm may accept. */ + private AllowedSolution allowed; + + /** + * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively) + */ + public BracketingNthOrderBrentSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double absoluteAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(absoluteAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolution.ANY_SIDE; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolution.ANY_SIDE; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolution.ANY_SIDE; + } + + /** Get the maximal order. + * @return maximal order + */ + public int getMaximalOrder() { + return maximalOrder; + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException, + NumberIsTooLargeException, + NoBracketingException { + // prepare arrays with the first points + final double[] x = new double[maximalOrder + 1]; + final double[] y = new double[maximalOrder + 1]; + x[0] = getMin(); + x[1] = getStartValue(); + x[2] = getMax(); + verifySequence(x[0], x[1], x[2]); + + // evaluate initial guess + y[1] = computeObjectiveValue(x[1]); + if (Precision.equals(y[1], 0.0, 1)) { + // return the initial guess if it is a perfect root. + return x[1]; + } + + // evaluate first endpoint + y[0] = computeObjectiveValue(x[0]); + if (Precision.equals(y[0], 0.0, 1)) { + // return the first endpoint if it is a perfect root. + return x[0]; + } + + int nbPoints; + int signChangeIndex; + if (y[0] * y[1] < 0) { + + // reduce interval if it brackets the root + nbPoints = 2; + signChangeIndex = 1; + + } else { + + // evaluate second endpoint + y[2] = computeObjectiveValue(x[2]); + if (Precision.equals(y[2], 0.0, 1)) { + // return the second endpoint if it is a perfect root. + return x[2]; + } + + if (y[1] * y[2] < 0) { + // use all computed point as a start sampling array for solving + nbPoints = 3; + signChangeIndex = 2; + } else { + throw new NoBracketingException(x[0], x[2], y[0], y[2]); + } + + } + + // prepare a work array for inverse polynomial interpolation + final double[] tmpX = new double[x.length]; + + // current tightest bracketing of the root + double xA = x[signChangeIndex - 1]; + double yA = y[signChangeIndex - 1]; + double absYA = FastMath.abs(yA); + int agingA = 0; + double xB = x[signChangeIndex]; + double yB = y[signChangeIndex]; + double absYB = FastMath.abs(yB); + int agingB = 0; + + // search loop + while (true) { + + // check convergence of bracketing interval + final double xTol = getAbsoluteAccuracy() + + getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB)); + if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) { + switch (allowed) { + case ANY_SIDE : + return absYA < absYB ? xA : xB; + case LEFT_SIDE : + return xA; + case RIGHT_SIDE : + return xB; + case BELOW_SIDE : + return (yA <= 0) ? xA : xB; + case ABOVE_SIDE : + return (yA < 0) ? xB : xA; + default : + // this should never happen + throw new MathInternalError(); + } + } + + // target for the next evaluation point + double targetY; + if (agingA >= MAXIMAL_AGING) { + // we keep updating the high bracket, try to compensate this + final int p = agingA - MAXIMAL_AGING; + final double weightA = (1 << p) - 1; + final double weightB = p + 1; + targetY = (weightA * yA - weightB * REDUCTION_FACTOR * yB) / (weightA + weightB); + } else if (agingB >= MAXIMAL_AGING) { + // we keep updating the low bracket, try to compensate this + final int p = agingB - MAXIMAL_AGING; + final double weightA = p + 1; + final double weightB = (1 << p) - 1; + targetY = (weightB * yB - weightA * REDUCTION_FACTOR * yA) / (weightA + weightB); + } else { + // bracketing is balanced, try to find the root itself + targetY = 0; + } + + // make a few attempts to guess a root, + double nextX; + int start = 0; + int end = nbPoints; + do { + + // guess a value for current target, using inverse polynomial interpolation + System.arraycopy(x, start, tmpX, start, end - start); + nextX = guessX(targetY, tmpX, y, start, end); + + if (!((nextX > xA) && (nextX < xB))) { + // the guessed root is not strictly inside of the tightest bracketing interval + + // the guessed root is either not strictly inside the interval or it + // is a NaN (which occurs when some sampling points share the same y) + // we try again with a lower interpolation order + if (signChangeIndex - start >= end - signChangeIndex) { + // we have more points before the sign change, drop the lowest point + ++start; + } else { + // we have more points after sign change, drop the highest point + --end; + } + + // we need to do one more attempt + nextX = Double.NaN; + + } + + } while (Double.isNaN(nextX) && (end - start > 1)); + + if (Double.isNaN(nextX)) { + // fall back to bisection + nextX = xA + 0.5 * (xB - xA); + start = signChangeIndex - 1; + end = signChangeIndex; + } + + // evaluate the function at the guessed root + final double nextY = computeObjectiveValue(nextX); + if (Precision.equals(nextY, 0.0, 1)) { + // we have found an exact root, since it is not an approximation + // we don't need to bother about the allowed solutions setting + return nextX; + } + + if ((nbPoints > 2) && (end - start != nbPoints)) { + + // we have been forced to ignore some points to keep bracketing, + // they are probably too far from the root, drop them from now on + nbPoints = end - start; + System.arraycopy(x, start, x, 0, nbPoints); + System.arraycopy(y, start, y, 0, nbPoints); + signChangeIndex -= start; + + } else if (nbPoints == x.length) { + + // we have to drop one point in order to insert the new one + nbPoints--; + + // keep the tightest bracketing interval as centered as possible + if (signChangeIndex >= (x.length + 1) / 2) { + // we drop the lowest point, we have to shift the arrays and the index + System.arraycopy(x, 1, x, 0, nbPoints); + System.arraycopy(y, 1, y, 0, nbPoints); + --signChangeIndex; + } + + } + + // insert the last computed point + //(by construction, we know it lies inside the tightest bracketing interval) + System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex); + x[signChangeIndex] = nextX; + System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex); + y[signChangeIndex] = nextY; + ++nbPoints; + + // update the bracketing interval + if (nextY * yA <= 0) { + // the sign change occurs before the inserted point + xB = nextX; + yB = nextY; + absYB = FastMath.abs(yB); + ++agingA; + agingB = 0; + } else { + // the sign change occurs after the inserted point + xA = nextX; + yA = nextY; + absYA = FastMath.abs(yA); + agingA = 0; + ++agingB; + + // update the sign change index + signChangeIndex++; + + } + + } + + } + + /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation. + * <p> + * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q + * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>), + * Q(y<sub>i</sub>) = x<sub>i</sub>. + * </p> + * @param targetY target value for y + * @param x reference points abscissas for interpolation, + * note that this array <em>is</em> modified during computation + * @param y reference points ordinates for interpolation + * @param start start index of the points to consider (inclusive) + * @param end end index of the points to consider (exclusive) + * @return guessed root (will be a NaN if two points share the same y) + */ + private double guessX(final double targetY, final double[] x, final double[] y, + final int start, final int end) { + + // compute Q Newton coefficients by divided differences + for (int i = start; i < end - 1; ++i) { + final int delta = i + 1 - start; + for (int j = end - 1; j > i; --j) { + x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]); + } + } + + // evaluate Q(targetY) + double x0 = 0; + for (int j = end - 1; j >= start; --j) { + x0 = x[j] + x0 * (targetY - y[j]); + } + + return x0; + + } + + /** {@inheritDoc} */ + public double solve(int maxEval, UnivariateFunction f, double min, + double max, AllowedSolution allowedSolution) + throws TooManyEvaluationsException, + NumberIsTooLargeException, + NoBracketingException { + this.allowed = allowedSolution; + return super.solve(maxEval, f, min, max); + } + + /** {@inheritDoc} */ + public double solve(int maxEval, UnivariateFunction f, double min, + double max, double startValue, + AllowedSolution allowedSolution) + throws TooManyEvaluationsException, + NumberIsTooLargeException, + NoBracketingException { + this.allowed = allowedSolution; + return super.solve(maxEval, f, min, max, startValue); + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java new file mode 100644 index 0000000..cf69410 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/BrentSolver.java @@ -0,0 +1,243 @@ +/* + * 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; +import org.apache.commons.math3.util.Precision; + +/** + * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html"> + * Brent algorithm</a> for finding zeros of real univariate functions. + * The function should be continuous but not necessarily smooth. + * The {@code solve} method returns a zero {@code x} of the function {@code f} + * in the given interval {@code [a, b]} to within a tolerance + * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and + * {@code t} is the absolute accuracy. + * <p>The given interval must bracket the root.</p> + * <p> + * The reference implementation is given in chapter 4 of + * <blockquote> + * <b>Algorithms for Minimization Without Derivatives</b>, + * <em>Richard P. Brent</em>, + * Dover, 2002 + * </blockquote> + * + * @see BaseAbstractUnivariateSolver + */ +public class BrentSolver extends AbstractUnivariateSolver { + + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver with default absolute accuracy (1e-6). + */ + public BrentSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public BrentSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public BrentSolver(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + * + * @see BaseAbstractUnivariateSolver#BaseAbstractUnivariateSolver(double,double,double) + */ + public BrentSolver(double relativeAccuracy, + double absoluteAccuracy, + double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws NoBracketingException, + TooManyEvaluationsException, + NumberIsTooLargeException { + double min = getMin(); + double max = getMax(); + final double initial = getStartValue(); + final double functionValueAccuracy = getFunctionValueAccuracy(); + + verifySequence(min, initial, max); + + // Return the initial guess if it is good enough. + double yInitial = computeObjectiveValue(initial); + if (FastMath.abs(yInitial) <= functionValueAccuracy) { + return initial; + } + + // Return the first endpoint if it is good enough. + double yMin = computeObjectiveValue(min); + if (FastMath.abs(yMin) <= functionValueAccuracy) { + return min; + } + + // Reduce interval if min and initial bracket the root. + if (yInitial * yMin < 0) { + return brent(min, initial, yMin, yInitial); + } + + // Return the second endpoint if it is good enough. + double yMax = computeObjectiveValue(max); + if (FastMath.abs(yMax) <= functionValueAccuracy) { + return max; + } + + // Reduce interval if initial and max bracket the root. + if (yInitial * yMax < 0) { + return brent(initial, max, yInitial, yMax); + } + + throw new NoBracketingException(min, max, yMin, yMax); + } + + /** + * Search for a zero inside the provided interval. + * This implementation is based on the algorithm described at page 58 of + * the book + * <blockquote> + * <b>Algorithms for Minimization Without Derivatives</b>, + * <it>Richard P. Brent</it>, + * Dover 0-486-41998-3 + * </blockquote> + * + * @param lo Lower bound of the search interval. + * @param hi Higher bound of the search interval. + * @param fLo Function value at the lower bound of the search interval. + * @param fHi Function value at the higher bound of the search interval. + * @return the value where the function is zero. + */ + private double brent(double lo, double hi, + double fLo, double fHi) { + double a = lo; + double fa = fLo; + double b = hi; + double fb = fHi; + double c = a; + double fc = fa; + double d = b - a; + double e = d; + + final double t = getAbsoluteAccuracy(); + final double eps = getRelativeAccuracy(); + + while (true) { + if (FastMath.abs(fc) < FastMath.abs(fb)) { + a = b; + b = c; + c = a; + fa = fb; + fb = fc; + fc = fa; + } + + final double tol = 2 * eps * FastMath.abs(b) + t; + final double m = 0.5 * (c - b); + + if (FastMath.abs(m) <= tol || + Precision.equals(fb, 0)) { + return b; + } + if (FastMath.abs(e) < tol || + FastMath.abs(fa) <= FastMath.abs(fb)) { + // Force bisection. + d = m; + e = d; + } else { + double s = fb / fa; + double p; + double q; + // The equality test (a == c) is intentional, + // it is part of the original Brent's method and + // it should NOT be replaced by proximity test. + if (a == c) { + // Linear interpolation. + p = 2 * m * s; + q = 1 - s; + } else { + // Inverse quadratic interpolation. + q = fa / fc; + final double r = fb / fc; + p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); + q = (q - 1) * (r - 1) * (s - 1); + } + if (p > 0) { + q = -q; + } else { + p = -p; + } + s = e; + e = d; + if (p >= 1.5 * m * q - FastMath.abs(tol * q) || + p >= FastMath.abs(0.5 * s * q)) { + // Inverse quadratic interpolation gives a value + // in the wrong direction, or progress is slow. + // Fall back to bisection. + d = m; + e = d; + } else { + d = p / q; + } + } + a = b; + fa = fb; + + if (FastMath.abs(d) > tol) { + b += d; + } else if (m > 0) { + b += tol; + } else { + b -= tol; + } + fb = computeObjectiveValue(b); + if ((fb > 0 && fc > 0) || + (fb <= 0 && fc <= 0)) { + c = a; + fc = fa; + d = b - a; + e = d; + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java new file mode 100644 index 0000000..b9ae158 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/DifferentiableUnivariateSolver.java @@ -0,0 +1,30 @@ +/* + * 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.analysis.DifferentiableUnivariateFunction; + + +/** + * Interface for (univariate real) rootfinding algorithms. + * Implementations will search for only one zero in the given interval. + * + * @deprecated as of 3.1, replaced by {@link UnivariateDifferentiableSolver} + */ +@Deprecated +public interface DifferentiableUnivariateSolver + extends BaseUnivariateSolver<DifferentiableUnivariateFunction> {} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/FieldBracketingNthOrderBrentSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/FieldBracketingNthOrderBrentSolver.java new file mode 100644 index 0000000..f0ca8b9 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/FieldBracketingNthOrderBrentSolver.java @@ -0,0 +1,446 @@ +/* + * 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.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.analysis.RealFieldUnivariateFunction; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.IntegerSequence; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; +import org.apache.commons.math3.util.Precision; + +/** + * This class implements a modification of the <a + * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>. + * <p> + * The changes with respect to the original Brent algorithm are: + * <ul> + * <li>the returned value is chosen in the current interval according + * to user specified {@link AllowedSolution}</li> + * <li>the maximal order for the invert polynomial root search is + * user-specified instead of being invert quadratic only</li> + * </ul><p> + * The given interval must bracket the root.</p> + * + * @param <T> the type of the field elements + * @since 3.6 + */ +public class FieldBracketingNthOrderBrentSolver<T extends RealFieldElement<T>> + implements BracketedRealFieldUnivariateSolver<T> { + + /** Maximal aging triggering an attempt to balance the bracketing interval. */ + private static final int MAXIMAL_AGING = 2; + + /** Field to which the elements belong. */ + private final Field<T> field; + + /** Maximal order. */ + private final int maximalOrder; + + /** Function value accuracy. */ + private final T functionValueAccuracy; + + /** Absolute accuracy. */ + private final T absoluteAccuracy; + + /** Relative accuracy. */ + private final T relativeAccuracy; + + /** Evaluations counter. */ + private IntegerSequence.Incrementor evaluations; + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public FieldBracketingNthOrderBrentSolver(final T relativeAccuracy, + final T absoluteAccuracy, + final T functionValueAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.field = relativeAccuracy.getField(); + this.maximalOrder = maximalOrder; + this.absoluteAccuracy = absoluteAccuracy; + this.relativeAccuracy = relativeAccuracy; + this.functionValueAccuracy = functionValueAccuracy; + this.evaluations = IntegerSequence.Incrementor.create(); + } + + /** Get the maximal order. + * @return maximal order + */ + public int getMaximalOrder() { + return maximalOrder; + } + + /** + * Get the maximal number of function evaluations. + * + * @return the maximal number of function evaluations. + */ + public int getMaxEvaluations() { + return evaluations.getMaximalCount(); + } + + /** + * Get the number of evaluations of the objective function. + * The number of evaluations corresponds to the last call to the + * {@code optimize} method. It is 0 if the method has not been + * called yet. + * + * @return the number of evaluations of the objective function. + */ + public int getEvaluations() { + return evaluations.getCount(); + } + + /** + * Get the absolute accuracy. + * @return absolute accuracy + */ + public T getAbsoluteAccuracy() { + return absoluteAccuracy; + } + + /** + * Get the relative accuracy. + * @return relative accuracy + */ + public T getRelativeAccuracy() { + return relativeAccuracy; + } + + /** + * Get the function accuracy. + * @return function accuracy + */ + public T getFunctionValueAccuracy() { + return functionValueAccuracy; + } + + /** + * Solve for a zero in the given interval. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return a value where the function is zero. + * @exception NullArgumentException if f is null. + * @exception NoBracketingException if root cannot be bracketed + */ + public T solve(final int maxEval, final RealFieldUnivariateFunction<T> f, + final T min, final T max, final AllowedSolution allowedSolution) + throws NullArgumentException, NoBracketingException { + return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution); + } + + /** + * Solve for a zero in the given interval, start at {@code startValue}. + * A solver may require that the interval brackets a single zero root. + * Solvers that do require bracketing should be able to handle the case + * where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @param allowedSolution The kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return a value where the function is zero. + * @exception NullArgumentException if f is null. + * @exception NoBracketingException if root cannot be bracketed + */ + public T solve(final int maxEval, final RealFieldUnivariateFunction<T> f, + final T min, final T max, final T startValue, + final AllowedSolution allowedSolution) + throws NullArgumentException, NoBracketingException { + + // Checks. + MathUtils.checkNotNull(f); + + // Reset. + evaluations = evaluations.withMaximalCount(maxEval).withStart(0); + T zero = field.getZero(); + T nan = zero.add(Double.NaN); + + // prepare arrays with the first points + final T[] x = MathArrays.buildArray(field, maximalOrder + 1); + final T[] y = MathArrays.buildArray(field, maximalOrder + 1); + x[0] = min; + x[1] = startValue; + x[2] = max; + + // evaluate initial guess + evaluations.increment(); + y[1] = f.value(x[1]); + if (Precision.equals(y[1].getReal(), 0.0, 1)) { + // return the initial guess if it is a perfect root. + return x[1]; + } + + // evaluate first endpoint + evaluations.increment(); + y[0] = f.value(x[0]); + if (Precision.equals(y[0].getReal(), 0.0, 1)) { + // return the first endpoint if it is a perfect root. + return x[0]; + } + + int nbPoints; + int signChangeIndex; + if (y[0].multiply(y[1]).getReal() < 0) { + + // reduce interval if it brackets the root + nbPoints = 2; + signChangeIndex = 1; + + } else { + + // evaluate second endpoint + evaluations.increment(); + y[2] = f.value(x[2]); + if (Precision.equals(y[2].getReal(), 0.0, 1)) { + // return the second endpoint if it is a perfect root. + return x[2]; + } + + if (y[1].multiply(y[2]).getReal() < 0) { + // use all computed point as a start sampling array for solving + nbPoints = 3; + signChangeIndex = 2; + } else { + throw new NoBracketingException(x[0].getReal(), x[2].getReal(), + y[0].getReal(), y[2].getReal()); + } + + } + + // prepare a work array for inverse polynomial interpolation + final T[] tmpX = MathArrays.buildArray(field, x.length); + + // current tightest bracketing of the root + T xA = x[signChangeIndex - 1]; + T yA = y[signChangeIndex - 1]; + T absXA = xA.abs(); + T absYA = yA.abs(); + int agingA = 0; + T xB = x[signChangeIndex]; + T yB = y[signChangeIndex]; + T absXB = xB.abs(); + T absYB = yB.abs(); + int agingB = 0; + + // search loop + while (true) { + + // check convergence of bracketing interval + T maxX = absXA.subtract(absXB).getReal() < 0 ? absXB : absXA; + T maxY = absYA.subtract(absYB).getReal() < 0 ? absYB : absYA; + final T xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX)); + if (xB.subtract(xA).subtract(xTol).getReal() <= 0 || + maxY.subtract(functionValueAccuracy).getReal() < 0) { + switch (allowedSolution) { + case ANY_SIDE : + return absYA.subtract(absYB).getReal() < 0 ? xA : xB; + case LEFT_SIDE : + return xA; + case RIGHT_SIDE : + return xB; + case BELOW_SIDE : + return yA.getReal() <= 0 ? xA : xB; + case ABOVE_SIDE : + return yA.getReal() < 0 ? xB : xA; + default : + // this should never happen + throw new MathInternalError(null); + } + } + + // target for the next evaluation point + T targetY; + if (agingA >= MAXIMAL_AGING) { + // we keep updating the high bracket, try to compensate this + targetY = yB.divide(16).negate(); + } else if (agingB >= MAXIMAL_AGING) { + // we keep updating the low bracket, try to compensate this + targetY = yA.divide(16).negate(); + } else { + // bracketing is balanced, try to find the root itself + targetY = zero; + } + + // make a few attempts to guess a root, + T nextX; + int start = 0; + int end = nbPoints; + do { + + // guess a value for current target, using inverse polynomial interpolation + System.arraycopy(x, start, tmpX, start, end - start); + nextX = guessX(targetY, tmpX, y, start, end); + + if (!((nextX.subtract(xA).getReal() > 0) && (nextX.subtract(xB).getReal() < 0))) { + // the guessed root is not strictly inside of the tightest bracketing interval + + // the guessed root is either not strictly inside the interval or it + // is a NaN (which occurs when some sampling points share the same y) + // we try again with a lower interpolation order + if (signChangeIndex - start >= end - signChangeIndex) { + // we have more points before the sign change, drop the lowest point + ++start; + } else { + // we have more points after sign change, drop the highest point + --end; + } + + // we need to do one more attempt + nextX = nan; + + } + + } while (Double.isNaN(nextX.getReal()) && (end - start > 1)); + + if (Double.isNaN(nextX.getReal())) { + // fall back to bisection + nextX = xA.add(xB.subtract(xA).divide(2)); + start = signChangeIndex - 1; + end = signChangeIndex; + } + + // evaluate the function at the guessed root + evaluations.increment(); + final T nextY = f.value(nextX); + if (Precision.equals(nextY.getReal(), 0.0, 1)) { + // we have found an exact root, since it is not an approximation + // we don't need to bother about the allowed solutions setting + return nextX; + } + + if ((nbPoints > 2) && (end - start != nbPoints)) { + + // we have been forced to ignore some points to keep bracketing, + // they are probably too far from the root, drop them from now on + nbPoints = end - start; + System.arraycopy(x, start, x, 0, nbPoints); + System.arraycopy(y, start, y, 0, nbPoints); + signChangeIndex -= start; + + } else if (nbPoints == x.length) { + + // we have to drop one point in order to insert the new one + nbPoints--; + + // keep the tightest bracketing interval as centered as possible + if (signChangeIndex >= (x.length + 1) / 2) { + // we drop the lowest point, we have to shift the arrays and the index + System.arraycopy(x, 1, x, 0, nbPoints); + System.arraycopy(y, 1, y, 0, nbPoints); + --signChangeIndex; + } + + } + + // insert the last computed point + //(by construction, we know it lies inside the tightest bracketing interval) + System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex); + x[signChangeIndex] = nextX; + System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex); + y[signChangeIndex] = nextY; + ++nbPoints; + + // update the bracketing interval + if (nextY.multiply(yA).getReal() <= 0) { + // the sign change occurs before the inserted point + xB = nextX; + yB = nextY; + absYB = yB.abs(); + ++agingA; + agingB = 0; + } else { + // the sign change occurs after the inserted point + xA = nextX; + yA = nextY; + absYA = yA.abs(); + agingA = 0; + ++agingB; + + // update the sign change index + signChangeIndex++; + + } + + } + + } + + /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation. + * <p> + * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q + * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>), + * Q(y<sub>i</sub>) = x<sub>i</sub>. + * </p> + * @param targetY target value for y + * @param x reference points abscissas for interpolation, + * note that this array <em>is</em> modified during computation + * @param y reference points ordinates for interpolation + * @param start start index of the points to consider (inclusive) + * @param end end index of the points to consider (exclusive) + * @return guessed root (will be a NaN if two points share the same y) + */ + private T guessX(final T targetY, final T[] x, final T[] y, + final int start, final int end) { + + // compute Q Newton coefficients by divided differences + for (int i = start; i < end - 1; ++i) { + final int delta = i + 1 - start; + for (int j = end - 1; j > i; --j) { + x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta])); + } + } + + // evaluate Q(targetY) + T x0 = field.getZero(); + for (int j = end - 1; j >= start; --j) { + x0 = x[j].add(x0.multiply(targetY.subtract(y[j]))); + } + + return x0; + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java new file mode 100644 index 0000000..bd3bc71 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/IllinoisSolver.java @@ -0,0 +1,82 @@ +/* + * 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; + + +/** + * Implements the <em>Illinois</em> method for root-finding (approximating + * a zero of a univariate real function). It is a modified + * {@link RegulaFalsiSolver <em>Regula Falsi</em>} method. + * + * <p>Like the <em>Regula Falsi</em> method, convergence is guaranteed by + * maintaining a bracketed solution. The <em>Illinois</em> method however, + * should converge much faster than the original <em>Regula Falsi</em> + * method. Furthermore, this implementation of the <em>Illinois</em> method + * should not suffer from the same implementation issues as the <em>Regula + * Falsi</em> method, which may fail to convergence in certain cases.</p> + * + * <p>The <em>Illinois</em> method assumes that the function is continuous, + * but not necessarily smooth.</p> + * + * <p>Implementation based on the following article: M. Dowell and P. Jarratt, + * <em>A modified regula falsi method for computing the root of an + * equation</em>, BIT Numerical Mathematics, volume 11, number 2, + * pages 168-174, Springer, 1971.</p> + * + * @since 3.0 + */ +public class IllinoisSolver extends BaseSecantSolver { + + /** Construct a solver with default accuracy (1e-6). */ + public IllinoisSolver() { + super(DEFAULT_ABSOLUTE_ACCURACY, Method.ILLINOIS); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public IllinoisSolver(final double absoluteAccuracy) { + super(absoluteAccuracy, Method.ILLINOIS); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public IllinoisSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy, Method.ILLINOIS); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Maximum function value error. + */ + public IllinoisSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.PEGASUS); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java new file mode 100644 index 0000000..5312dac --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/LaguerreSolver.java @@ -0,0 +1,440 @@ +/* + * 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.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.complex.Complex; +import org.apache.commons.math3.complex.ComplexUtils; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.FastMath; + +/** + * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html"> + * Laguerre's Method</a> for root finding of real coefficient polynomials. + * For reference, see + * <blockquote> + * <b>A First Course in Numerical Analysis</b>, + * ISBN 048641454X, chapter 8. + * </blockquote> + * Laguerre's method is global in the sense that it can start with any initial + * approximation and be able to solve all roots from that point. + * The algorithm requires a bracketing condition. + * + * @since 1.2 + */ +public class LaguerreSolver extends AbstractPolynomialSolver { + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + /** Complex solver. */ + private final ComplexSolver complexSolver = new ComplexSolver(); + + /** + * Construct a solver with default accuracy (1e-6). + */ + public LaguerreSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public LaguerreSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public LaguerreSolver(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + */ + public LaguerreSolver(double relativeAccuracy, + double absoluteAccuracy, + double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + public 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); + + // Return the initial guess if it is good enough. + final double yInitial = computeObjectiveValue(initial); + if (FastMath.abs(yInitial) <= functionValueAccuracy) { + return initial; + } + + // Return the first endpoint if it is good enough. + final double yMin = computeObjectiveValue(min); + if (FastMath.abs(yMin) <= functionValueAccuracy) { + return min; + } + + // Reduce interval if min and initial bracket the root. + if (yInitial * yMin < 0) { + return laguerre(min, initial, yMin, yInitial); + } + + // Return the second endpoint if it is good enough. + final double yMax = computeObjectiveValue(max); + if (FastMath.abs(yMax) <= functionValueAccuracy) { + return max; + } + + // Reduce interval if initial and max bracket the root. + if (yInitial * yMax < 0) { + return laguerre(initial, max, yInitial, yMax); + } + + throw new NoBracketingException(min, max, yMin, yMax); + } + + /** + * Find a real root in the given interval. + * + * Despite the bracketing condition, the root returned by + * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may + * not be a real zero inside {@code [min, max]}. + * For example, <code> p(x) = x<sup>3</sup> + 1, </code> + * with {@code min = -2}, {@code max = 2}, {@code initial = 0}. + * When it occurs, this code calls + * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)} + * in order to obtain all roots and picks up one real root. + * + * @param lo Lower bound of the search interval. + * @param hi Higher bound of the search interval. + * @param fLo Function value at the lower bound of the search interval. + * @param fHi Function value at the higher bound of the search interval. + * @return the point at which the function value is zero. + * @deprecated This method should not be part of the public API: It will + * be made private in version 4.0. + */ + @Deprecated + public double laguerre(double lo, double hi, + double fLo, double fHi) { + final Complex c[] = ComplexUtils.convertToComplex(getCoefficients()); + + final Complex initial = new Complex(0.5 * (lo + hi), 0); + final Complex z = complexSolver.solve(c, initial); + if (complexSolver.isRoot(lo, hi, z)) { + return z.getReal(); + } else { + double r = Double.NaN; + // Solve all roots and select the one we are seeking. + Complex[] root = complexSolver.solveAll(c, initial); + for (int i = 0; i < root.length; i++) { + if (complexSolver.isRoot(lo, hi, root[i])) { + r = root[i].getReal(); + break; + } + } + return r; + } + } + + /** + * Find all complex roots for the polynomial with the given + * coefficients, starting from the given initial value. + * <p> + * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> + * + * @param coefficients Polynomial coefficients. + * @param initial Start value. + * @return the full set of complex roots of the polynomial + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded when solving for one of the roots + * @throws NullArgumentException if the {@code coefficients} is + * {@code null}. + * @throws NoDataException if the {@code coefficients} array is empty. + * @since 3.1 + */ + public Complex[] solveAllComplex(double[] coefficients, + double initial) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + return solveAllComplex(coefficients, initial, Integer.MAX_VALUE); + } + + /** + * Find all complex roots for the polynomial with the given + * coefficients, starting from the given initial value. + * <p> + * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> + * + * @param coefficients polynomial coefficients + * @param initial start value + * @param maxEval maximum number of evaluations + * @return the full set of complex roots of the polynomial + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded when solving for one of the roots + * @throws NullArgumentException if the {@code coefficients} is + * {@code null} + * @throws NoDataException if the {@code coefficients} array is empty + * @since 3.5 + */ + public Complex[] solveAllComplex(double[] coefficients, + double initial, int maxEval) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + setup(maxEval, + new PolynomialFunction(coefficients), + Double.NEGATIVE_INFINITY, + Double.POSITIVE_INFINITY, + initial); + return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients), + new Complex(initial, 0d)); + } + + /** + * Find a complex root for the polynomial with the given coefficients, + * starting from the given initial value. + * <p> + * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> + * + * @param coefficients Polynomial coefficients. + * @param initial Start value. + * @return a complex root of the polynomial + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded. + * @throws NullArgumentException if the {@code coefficients} is + * {@code null}. + * @throws NoDataException if the {@code coefficients} array is empty. + * @since 3.1 + */ + public Complex solveComplex(double[] coefficients, + double initial) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + return solveComplex(coefficients, initial, Integer.MAX_VALUE); + } + + /** + * Find a complex root for the polynomial with the given coefficients, + * starting from the given initial value. + * <p> + * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p> + * + * @param coefficients polynomial coefficients + * @param initial start value + * @param maxEval maximum number of evaluations + * @return a complex root of the polynomial + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded + * @throws NullArgumentException if the {@code coefficients} is + * {@code null} + * @throws NoDataException if the {@code coefficients} array is empty + * @since 3.1 + */ + public Complex solveComplex(double[] coefficients, + double initial, int maxEval) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + setup(maxEval, + new PolynomialFunction(coefficients), + Double.NEGATIVE_INFINITY, + Double.POSITIVE_INFINITY, + initial); + return complexSolver.solve(ComplexUtils.convertToComplex(coefficients), + new Complex(initial, 0d)); + } + + /** + * Class for searching all (complex) roots. + */ + private class ComplexSolver { + /** + * Check whether the given complex root is actually a real zero + * in the given interval, within the solver tolerance level. + * + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param z Complex root. + * @return {@code true} if z is a real zero. + */ + public boolean isRoot(double min, double max, Complex z) { + if (isSequence(min, z.getReal(), max)) { + double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy()); + return (FastMath.abs(z.getImaginary()) <= tolerance) || + (z.abs() <= getFunctionValueAccuracy()); + } + return false; + } + + /** + * Find all complex roots for the polynomial with the given + * coefficients, starting from the given initial value. + * + * @param coefficients Polynomial coefficients. + * @param initial Start value. + * @return the point at which the function value is zero. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded. + * @throws NullArgumentException if the {@code coefficients} is + * {@code null}. + * @throws NoDataException if the {@code coefficients} array is empty. + */ + public Complex[] solveAll(Complex coefficients[], Complex initial) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + if (coefficients == null) { + throw new NullArgumentException(); + } + final int n = coefficients.length - 1; + if (n == 0) { + throw new NoDataException(LocalizedFormats.POLYNOMIAL); + } + // Coefficients for deflated polynomial. + final Complex c[] = new Complex[n + 1]; + for (int i = 0; i <= n; i++) { + c[i] = coefficients[i]; + } + + // Solve individual roots successively. + final Complex root[] = new Complex[n]; + for (int i = 0; i < n; i++) { + final Complex subarray[] = new Complex[n - i + 1]; + System.arraycopy(c, 0, subarray, 0, subarray.length); + root[i] = solve(subarray, initial); + // Polynomial deflation using synthetic division. + Complex newc = c[n - i]; + Complex oldc = null; + for (int j = n - i - 1; j >= 0; j--) { + oldc = c[j]; + c[j] = newc; + newc = oldc.add(newc.multiply(root[i])); + } + } + + return root; + } + + /** + * Find a complex root for the polynomial with the given coefficients, + * starting from the given initial value. + * + * @param coefficients Polynomial coefficients. + * @param initial Start value. + * @return the point at which the function value is zero. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum number of evaluations is exceeded. + * @throws NullArgumentException if the {@code coefficients} is + * {@code null}. + * @throws NoDataException if the {@code coefficients} array is empty. + */ + public Complex solve(Complex coefficients[], Complex initial) + throws NullArgumentException, + NoDataException, + TooManyEvaluationsException { + if (coefficients == null) { + throw new NullArgumentException(); + } + + final int n = coefficients.length - 1; + if (n == 0) { + throw new NoDataException(LocalizedFormats.POLYNOMIAL); + } + + final double absoluteAccuracy = getAbsoluteAccuracy(); + final double relativeAccuracy = getRelativeAccuracy(); + final double functionValueAccuracy = getFunctionValueAccuracy(); + + final Complex nC = new Complex(n, 0); + final Complex n1C = new Complex(n - 1, 0); + + Complex z = initial; + Complex oldz = new Complex(Double.POSITIVE_INFINITY, + Double.POSITIVE_INFINITY); + while (true) { + // Compute pv (polynomial value), dv (derivative value), and + // d2v (second derivative value) simultaneously. + Complex pv = coefficients[n]; + Complex dv = Complex.ZERO; + Complex d2v = Complex.ZERO; + for (int j = n-1; j >= 0; j--) { + d2v = dv.add(z.multiply(d2v)); + dv = pv.add(z.multiply(dv)); + pv = coefficients[j].add(z.multiply(pv)); + } + d2v = d2v.multiply(new Complex(2.0, 0.0)); + + // Check for convergence. + final double tolerance = FastMath.max(relativeAccuracy * z.abs(), + absoluteAccuracy); + if ((z.subtract(oldz)).abs() <= tolerance) { + return z; + } + if (pv.abs() <= functionValueAccuracy) { + return z; + } + + // Now pv != 0, calculate the new approximation. + final Complex G = dv.divide(pv); + final Complex G2 = G.multiply(G); + final Complex H = G2.subtract(d2v.divide(pv)); + final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2)); + // Choose a denominator larger in magnitude. + final Complex deltaSqrt = delta.sqrt(); + final Complex dplus = G.add(deltaSqrt); + final Complex dminus = G.subtract(deltaSqrt); + final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus; + // Perturb z if denominator is zero, for instance, + // p(x) = x^3 + 1, z = 0. + if (denominator.equals(new Complex(0.0, 0.0))) { + z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy)); + oldz = new Complex(Double.POSITIVE_INFINITY, + Double.POSITIVE_INFINITY); + } else { + oldz = z; + z = z.subtract(nC.divide(denominator)); + } + incrementEvaluationCount(); + } + } + } +} 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; + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java new file mode 100644 index 0000000..864cfd5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/MullerSolver2.java @@ -0,0 +1,168 @@ +/* + * 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> + * Except for the initial [min, max], it does not require bracketing + * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If a complex + * number arises in the computation, we simply use its modulus as a real + * approximation.</p> + * <p> + * Because the interval may not be bracketing, the bisection alternative is + * not applicable here. However in practice our treatment usually works + * well, especially near real zeroes where the imaginary part of the complex + * approximation is often negligible.</p> + * <p> + * The formulas here do not use divided differences directly.</p> + * + * @since 1.2 + * @see MullerSolver + */ +public class MullerSolver2 extends AbstractUnivariateSolver { + + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver with default accuracy (1e-6). + */ + public MullerSolver2() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public MullerSolver2(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public MullerSolver2(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException, + NumberIsTooLargeException, + NoBracketingException { + final double min = getMin(); + final double max = getMax(); + + verifyInterval(min, max); + + final double relativeAccuracy = getRelativeAccuracy(); + final double absoluteAccuracy = getAbsoluteAccuracy(); + final double functionValueAccuracy = getFunctionValueAccuracy(); + + // x2 is the last root approximation + // x is the new approximation and new x2 for next round + // x0 < x1 < x2 does not hold here + + double x0 = min; + double y0 = computeObjectiveValue(x0); + if (FastMath.abs(y0) < functionValueAccuracy) { + return x0; + } + double x1 = max; + double y1 = computeObjectiveValue(x1); + if (FastMath.abs(y1) < functionValueAccuracy) { + return x1; + } + + if(y0 * y1 > 0) { + throw new NoBracketingException(x0, x1, y0, y1); + } + + double x2 = 0.5 * (x0 + x1); + double y2 = computeObjectiveValue(x2); + + double oldx = Double.POSITIVE_INFINITY; + while (true) { + // quadratic interpolation through x0, x1, x2 + final double q = (x2 - x1) / (x1 - x0); + final double a = q * (y2 - (1 + q) * y1 + q * y0); + final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; + final double c = (1 + q) * y2; + final double delta = b * b - 4 * a * c; + double x; + final double denominator; + if (delta >= 0.0) { + // choose a denominator larger in magnitude + double dplus = b + FastMath.sqrt(delta); + double dminus = b - FastMath.sqrt(delta); + denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus; + } else { + // take the modulus of (B +/- FastMath.sqrt(delta)) + denominator = FastMath.sqrt(b * b - delta); + } + if (denominator != 0) { + x = x2 - 2.0 * c * (x2 - x1) / denominator; + // perturb x if it exactly coincides with x1 or x2 + // the equality tests here are intentional + while (x == x1 || x == x2) { + x += absoluteAccuracy; + } + } else { + // extremely rare case, get a random number to skip it + x = min + FastMath.random() * (max - min); + oldx = Double.POSITIVE_INFINITY; + } + 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; + } + + // prepare the next iteration + x0 = x1; + y0 = y1; + x1 = x2; + y1 = y2; + x2 = x; + y2 = y; + oldx = x; + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java new file mode 100644 index 0000000..4cf2688 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonRaphsonSolver.java @@ -0,0 +1,92 @@ +/* + * 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.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html"> + * Newton's Method</a> for finding zeros of real univariate differentiable + * functions. + * + * @since 3.1 + */ +public class NewtonRaphsonSolver extends AbstractUnivariateDifferentiableSolver { + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver. + */ + public NewtonRaphsonSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public NewtonRaphsonSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + + /** + * Find a zero near the midpoint of {@code min} and {@code max}. + * + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param maxEval Maximum number of evaluations. + * @return the value where the function is zero. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum evaluation count is exceeded. + * @throws org.apache.commons.math3.exception.NumberIsTooLargeException + * if {@code min >= max}. + */ + @Override + public double solve(int maxEval, final UnivariateDifferentiableFunction f, + final double min, final double max) + throws TooManyEvaluationsException { + return super.solve(maxEval, f, UnivariateSolverUtils.midpoint(min, max)); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException { + final double startValue = getStartValue(); + final double absoluteAccuracy = getAbsoluteAccuracy(); + + double x0 = startValue; + double x1; + while (true) { + final DerivativeStructure y0 = computeObjectiveValueAndDerivative(x0); + x1 = x0 - (y0.getValue() / y0.getPartialDerivative(1)); + if (FastMath.abs(x1 - x0) <= absoluteAccuracy) { + return x1; + } + + x0 = x1; + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java new file mode 100644 index 0000000..3ba7bf2 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/NewtonSolver.java @@ -0,0 +1,92 @@ +/* + * 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.analysis.DifferentiableUnivariateFunction; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html"> + * Newton's Method</a> for finding zeros of real univariate functions. + * <p> + * The function should be continuous but not necessarily smooth.</p> + * + * @deprecated as of 3.1, replaced by {@link NewtonRaphsonSolver} + */ +@Deprecated +public class NewtonSolver extends AbstractDifferentiableUnivariateSolver { + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver. + */ + public NewtonSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public NewtonSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + + /** + * Find a zero near the midpoint of {@code min} and {@code max}. + * + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param maxEval Maximum number of evaluations. + * @return the value where the function is zero. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximum evaluation count is exceeded. + * @throws org.apache.commons.math3.exception.NumberIsTooLargeException + * if {@code min >= max}. + */ + @Override + public double solve(int maxEval, final DifferentiableUnivariateFunction f, + final double min, final double max) + throws TooManyEvaluationsException { + return super.solve(maxEval, f, UnivariateSolverUtils.midpoint(min, max)); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException { + final double startValue = getStartValue(); + final double absoluteAccuracy = getAbsoluteAccuracy(); + + double x0 = startValue; + double x1; + while (true) { + x1 = x0 - (computeObjectiveValue(x0) / computeDerivativeObjectiveValue(x0)); + if (FastMath.abs(x1 - x0) <= absoluteAccuracy) { + return x1; + } + + x0 = x1; + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java new file mode 100644 index 0000000..0d80895 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/PegasusSolver.java @@ -0,0 +1,84 @@ +/* + * 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; + +/** + * Implements the <em>Pegasus</em> method for root-finding (approximating + * a zero of a univariate real function). It is a modified + * {@link RegulaFalsiSolver <em>Regula Falsi</em>} method. + * + * <p>Like the <em>Regula Falsi</em> method, convergence is guaranteed by + * maintaining a bracketed solution. The <em>Pegasus</em> method however, + * should converge much faster than the original <em>Regula Falsi</em> + * method. Furthermore, this implementation of the <em>Pegasus</em> method + * should not suffer from the same implementation issues as the <em>Regula + * Falsi</em> method, which may fail to convergence in certain cases. Also, + * the <em>Pegasus</em> method should converge faster than the + * {@link IllinoisSolver <em>Illinois</em>} method, another <em>Regula + * Falsi</em>-based method.</p> + * + * <p>The <em>Pegasus</em> method assumes that the function is continuous, + * but not necessarily smooth.</p> + * + * <p>Implementation based on the following article: M. Dowell and P. Jarratt, + * <em>The "Pegasus" method for computing the root of an equation</em>, + * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer, + * 1972.</p> + * + * @since 3.0 + */ +public class PegasusSolver extends BaseSecantSolver { + + /** Construct a solver with default accuracy (1e-6). */ + public PegasusSolver() { + super(DEFAULT_ABSOLUTE_ACCURACY, Method.PEGASUS); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public PegasusSolver(final double absoluteAccuracy) { + super(absoluteAccuracy, Method.PEGASUS); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public PegasusSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy, Method.PEGASUS); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Maximum function value error. + */ + public PegasusSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.PEGASUS); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java new file mode 100644 index 0000000..c21f076 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/PolynomialSolver.java @@ -0,0 +1,28 @@ +/* + * 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.analysis.polynomials.PolynomialFunction; + +/** + * Interface for (polynomial) root-finding algorithms. + * Implementations will search for only one zero in the given interval. + * + * @since 3.0 + */ +public interface PolynomialSolver + extends BaseUnivariateSolver<PolynomialFunction> {} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java new file mode 100644 index 0000000..cfb7055 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/RegulaFalsiSolver.java @@ -0,0 +1,94 @@ +/* + * 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; + +/** + * Implements the <em>Regula Falsi</em> or <em>False position</em> method for + * root-finding (approximating a zero of a univariate real function). It is a + * modified {@link SecantSolver <em>Secant</em>} method. + * + * <p>The <em>Regula Falsi</em> method is included for completeness, for + * testing purposes, for educational purposes, for comparison to other + * algorithms, etc. It is however <strong>not</strong> intended to be used + * for actual problems, as one of the bounds often remains fixed, resulting + * in very slow convergence. Instead, one of the well-known modified + * <em>Regula Falsi</em> algorithms can be used ({@link IllinoisSolver + * <em>Illinois</em>} or {@link PegasusSolver <em>Pegasus</em>}). These two + * algorithms solve the fundamental issues of the original <em>Regula + * Falsi</em> algorithm, and greatly out-performs it for most, if not all, + * (practical) functions. + * + * <p>Unlike the <em>Secant</em> method, the <em>Regula Falsi</em> guarantees + * convergence, by maintaining a bracketed solution. Note however, that due to + * the finite/limited precision of Java's {@link Double double} type, which is + * used in this implementation, the algorithm may get stuck in a situation + * where it no longer makes any progress. Such cases are detected and result + * in a {@code ConvergenceException} exception being thrown. In other words, + * the algorithm theoretically guarantees convergence, but the implementation + * does not.</p> + * + * <p>The <em>Regula Falsi</em> method assumes that the function is continuous, + * but not necessarily smooth.</p> + * + * <p>Implementation based on the following article: M. Dowell and P. Jarratt, + * <em>A modified regula falsi method for computing the root of an + * equation</em>, BIT Numerical Mathematics, volume 11, number 2, + * pages 168-174, Springer, 1971.</p> + * + * @since 3.0 + */ +public class RegulaFalsiSolver extends BaseSecantSolver { + + /** Construct a solver with default accuracy (1e-6). */ + public RegulaFalsiSolver() { + super(DEFAULT_ABSOLUTE_ACCURACY, Method.REGULA_FALSI); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public RegulaFalsiSolver(final double absoluteAccuracy) { + super(absoluteAccuracy, Method.REGULA_FALSI); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public RegulaFalsiSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy, Method.REGULA_FALSI); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Maximum function value error. + */ + public RegulaFalsiSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy) { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, Method.REGULA_FALSI); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java new file mode 100644 index 0000000..d83f595 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/RiddersSolver.java @@ -0,0 +1,142 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html"> + * Ridders' Method</a> for root finding of real univariate functions. For + * reference, see C. Ridders, <i>A new algorithm for computing a single root + * of a real continuous function </i>, IEEE Transactions on Circuits and + * Systems, 26 (1979), 979 - 980. + * <p> + * The function should be continuous but not necessarily smooth.</p> + * + * @since 1.2 + */ +public class RiddersSolver extends AbstractUnivariateSolver { + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** + * Construct a solver with default accuracy (1e-6). + */ + public RiddersSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY); + } + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + */ + public RiddersSolver(double absoluteAccuracy) { + super(absoluteAccuracy); + } + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + */ + public RiddersSolver(double relativeAccuracy, + double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() + throws TooManyEvaluationsException, + NoBracketingException { + double min = getMin(); + double max = getMax(); + // [x1, x2] is the bracketing interval in each iteration + // x3 is the midpoint of [x1, x2] + // x is the new root approximation and an endpoint of the new interval + double x1 = min; + double y1 = computeObjectiveValue(x1); + double x2 = max; + double y2 = computeObjectiveValue(x2); + + // check for zeros before verifying bracketing + if (y1 == 0) { + return min; + } + if (y2 == 0) { + return max; + } + verifyBracketing(min, max); + + final double absoluteAccuracy = getAbsoluteAccuracy(); + final double functionValueAccuracy = getFunctionValueAccuracy(); + final double relativeAccuracy = getRelativeAccuracy(); + + double oldx = Double.POSITIVE_INFINITY; + while (true) { + // calculate the new root approximation + final double x3 = 0.5 * (x1 + x2); + final double y3 = computeObjectiveValue(x3); + if (FastMath.abs(y3) <= functionValueAccuracy) { + return x3; + } + final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing + final double correction = (FastMath.signum(y2) * FastMath.signum(y3)) * + (x3 - x1) / FastMath.sqrt(delta); + final double x = x3 - correction; // correction != 0 + final double y = computeObjectiveValue(x); + + // check for convergence + final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy); + if (FastMath.abs(x - oldx) <= tolerance) { + return x; + } + if (FastMath.abs(y) <= functionValueAccuracy) { + return x; + } + + // prepare the new interval for next iteration + // Ridders' method guarantees x1 < x < x2 + if (correction > 0.0) { // x1 < x < x3 + if (FastMath.signum(y1) + FastMath.signum(y) == 0.0) { + x2 = x; + y2 = y; + } else { + x1 = x; + x2 = x3; + y1 = y; + y2 = y3; + } + } else { // x3 < x < x2 + if (FastMath.signum(y2) + FastMath.signum(y) == 0.0) { + x1 = x; + y1 = y; + } else { + x1 = x3; + x2 = x; + y1 = y3; + y2 = y; + } + } + oldx = x; + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/SecantSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/SecantSolver.java new file mode 100644 index 0000000..d866cf8 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/SecantSolver.java @@ -0,0 +1,135 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Implements the <em>Secant</em> method for root-finding (approximating a + * zero of a univariate real function). The solution that is maintained is + * not bracketed, and as such convergence is not guaranteed. + * + * <p>Implementation based on the following article: M. Dowell and P. Jarratt, + * <em>A modified regula falsi method for computing the root of an + * equation</em>, BIT Numerical Mathematics, volume 11, number 2, + * pages 168-174, Springer, 1971.</p> + * + * <p>Note that since release 3.0 this class implements the actual + * <em>Secant</em> algorithm, and not a modified one. As such, the 3.0 version + * is not backwards compatible with previous versions. To use an algorithm + * similar to the pre-3.0 releases, use the + * {@link IllinoisSolver <em>Illinois</em>} algorithm or the + * {@link PegasusSolver <em>Pegasus</em>} algorithm.</p> + * + */ +public class SecantSolver extends AbstractUnivariateSolver { + + /** Default absolute accuracy. */ + protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** Construct a solver with default accuracy (1e-6). */ + public SecantSolver() { + super(DEFAULT_ABSOLUTE_ACCURACY); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy absolute accuracy + */ + public SecantSolver(final double absoluteAccuracy) { + super(absoluteAccuracy); + } + + /** + * Construct a solver. + * + * @param relativeAccuracy relative accuracy + * @param absoluteAccuracy absolute accuracy + */ + public SecantSolver(final double relativeAccuracy, + final double absoluteAccuracy) { + super(relativeAccuracy, absoluteAccuracy); + } + + /** {@inheritDoc} */ + @Override + protected final double doSolve() + throws TooManyEvaluationsException, + NoBracketingException { + // Get initial solution + double x0 = getMin(); + double x1 = getMax(); + double f0 = computeObjectiveValue(x0); + double f1 = computeObjectiveValue(x1); + + // If one of the bounds is the exact root, return it. Since these are + // not under-approximations or over-approximations, we can return them + // regardless of the allowed solutions. + if (f0 == 0.0) { + return x0; + } + if (f1 == 0.0) { + return x1; + } + + // Verify bracketing of initial solution. + verifyBracketing(x0, x1); + + // Get accuracies. + final double ftol = getFunctionValueAccuracy(); + final double atol = getAbsoluteAccuracy(); + final double rtol = getRelativeAccuracy(); + + // Keep finding better approximations. + while (true) { + // Calculate the next approximation. + 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, + // we can return it regardless of the allowed solutions. + if (fx == 0.0) { + return x; + } + + // Update the bounds with the new approximation. + x0 = x1; + f0 = f1; + x1 = x; + f1 = fx; + + // If the function value of the last approximation is too small, + // given the function value accuracy, then we can't get closer to + // the root than we already are. + if (FastMath.abs(f1) <= ftol) { + return x1; + } + + // If the current interval is within the given accuracies, we + // are satisfied with the current approximation. + if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol)) { + return x1; + } + } + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateDifferentiableSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateDifferentiableSolver.java new file mode 100644 index 0000000..82bbead --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateDifferentiableSolver.java @@ -0,0 +1,29 @@ +/* + * 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.analysis.differentiation.UnivariateDifferentiableFunction; + + +/** + * Interface for (univariate real) rootfinding algorithms. + * Implementations will search for only one zero in the given interval. + * + * @since 3.1 + */ +public interface UnivariateDifferentiableSolver + extends BaseUnivariateSolver<UnivariateDifferentiableFunction> {} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolver.java b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolver.java new file mode 100644 index 0000000..484e67a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolver.java @@ -0,0 +1,28 @@ +/* + * 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.analysis.UnivariateFunction; + + +/** + * Interface for (univariate real) root-finding algorithms. + * Implementations will search for only one zero in the given interval. + * + */ +public interface UnivariateSolver + extends BaseUnivariateSolver<UnivariateFunction> {} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java new file mode 100644 index 0000000..71b34c5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java @@ -0,0 +1,467 @@ +/* + * 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.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.FastMath; + +/** + * Utility routines for {@link UnivariateSolver} objects. + * + */ +public class UnivariateSolverUtils { + /** + * Class contains only static methods. + */ + private UnivariateSolverUtils() {} + + /** + * Convenience method to find a zero of a univariate real function. A default + * solver is used. + * + * @param function Function. + * @param x0 Lower bound for the interval. + * @param x1 Upper bound for the interval. + * @return a value where the function is zero. + * @throws NoBracketingException if the function has the same sign at the + * endpoints. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static double solve(UnivariateFunction function, double x0, double x1) + throws NullArgumentException, + NoBracketingException { + if (function == null) { + throw new NullArgumentException(LocalizedFormats.FUNCTION); + } + final UnivariateSolver solver = new BrentSolver(); + return solver.solve(Integer.MAX_VALUE, function, x0, x1); + } + + /** + * Convenience method to find a zero of a univariate real function. A default + * solver is used. + * + * @param function Function. + * @param x0 Lower bound for the interval. + * @param x1 Upper bound for the interval. + * @param absoluteAccuracy Accuracy to be used by the solver. + * @return a value where the function is zero. + * @throws NoBracketingException if the function has the same sign at the + * endpoints. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static double solve(UnivariateFunction function, + double x0, double x1, + double absoluteAccuracy) + throws NullArgumentException, + NoBracketingException { + if (function == null) { + throw new NullArgumentException(LocalizedFormats.FUNCTION); + } + final UnivariateSolver solver = new BrentSolver(absoluteAccuracy); + return solver.solve(Integer.MAX_VALUE, function, x0, x1); + } + + /** + * Force a root found by a non-bracketing solver to lie on a specified side, + * as if the solver were a bracketing one. + * + * @param maxEval maximal number of new evaluations of the function + * (evaluations already done for finding the root should have already been subtracted + * from this number) + * @param f function to solve + * @param bracketing bracketing solver to use for shifting the root + * @param baseRoot original root found by a previous non-bracketing solver + * @param min minimal bound of the search interval + * @param max maximal bound of the search interval + * @param allowedSolution the kind of solutions that the root-finding algorithm may + * accept as solutions. + * @return a root approximation, on the specified side of the exact root + * @throws NoBracketingException if the function has the same sign at the + * endpoints. + */ + public static double forceSide(final int maxEval, final UnivariateFunction f, + final BracketedUnivariateSolver<UnivariateFunction> bracketing, + final double baseRoot, final double min, final double max, + final AllowedSolution allowedSolution) + throws NoBracketingException { + + if (allowedSolution == AllowedSolution.ANY_SIDE) { + // no further bracketing required + return baseRoot; + } + + // find a very small interval bracketing the root + final double step = FastMath.max(bracketing.getAbsoluteAccuracy(), + FastMath.abs(baseRoot * bracketing.getRelativeAccuracy())); + double xLo = FastMath.max(min, baseRoot - step); + double fLo = f.value(xLo); + double xHi = FastMath.min(max, baseRoot + step); + double fHi = f.value(xHi); + int remainingEval = maxEval - 2; + while (remainingEval > 0) { + + if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) { + // compute the root on the selected side + return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution); + } + + // try increasing the interval + boolean changeLo = false; + boolean changeHi = false; + if (fLo < fHi) { + // increasing function + if (fLo >= 0) { + changeLo = true; + } else { + changeHi = true; + } + } else if (fLo > fHi) { + // decreasing function + if (fLo <= 0) { + changeLo = true; + } else { + changeHi = true; + } + } else { + // unknown variation + changeLo = true; + changeHi = true; + } + + // update the lower bound + if (changeLo) { + xLo = FastMath.max(min, xLo - step); + fLo = f.value(xLo); + remainingEval--; + } + + // update the higher bound + if (changeHi) { + xHi = FastMath.min(max, xHi + step); + fHi = f.value(xHi); + remainingEval--; + } + + } + + throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING, + xLo, xHi, fLo, fHi, + maxEval - remainingEval, maxEval, baseRoot, + min, max); + + } + + /** + * This method simply calls {@link #bracket(UnivariateFunction, double, double, double, + * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)} + * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}. + * <p> + * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE} + * iterations to throw a {@code ConvergenceException.} Unless you are + * confident that there is a root between {@code lowerBound} and + * {@code upperBound} near {@code initial}, it is better to use + * {@link #bracket(UnivariateFunction, double, double, double, double,double, int) + * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}, + * explicitly specifying the maximum number of iterations.</p> + * + * @param function Function. + * @param initial Initial midpoint of interval being expanded to + * bracket a root. + * @param lowerBound Lower bound (a is never lower than this value) + * @param upperBound Upper bound (b never is greater than this + * value). + * @return a two-element array holding a and b. + * @throws NoBracketingException if a root cannot be bracketted. + * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static double[] bracket(UnivariateFunction function, + double initial, + double lowerBound, double upperBound) + throws NullArgumentException, + NotStrictlyPositiveException, + NoBracketingException { + return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, Integer.MAX_VALUE); + } + + /** + * This method simply calls {@link #bracket(UnivariateFunction, double, double, double, + * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)} + * with {@code q} and {@code r} set to 1.0. + * @param function Function. + * @param initial Initial midpoint of interval being expanded to + * bracket a root. + * @param lowerBound Lower bound (a is never lower than this value). + * @param upperBound Upper bound (b never is greater than this + * value). + * @param maximumIterations Maximum number of iterations to perform + * @return a two element array holding a and b. + * @throws NoBracketingException if the algorithm fails to find a and b + * satisfying the desired conditions. + * @throws NotStrictlyPositiveException if {@code maximumIterations <= 0}. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static double[] bracket(UnivariateFunction function, + double initial, + double lowerBound, double upperBound, + int maximumIterations) + throws NullArgumentException, + NotStrictlyPositiveException, + NoBracketingException { + return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations); + } + + /** + * This method attempts to find two values a and b satisfying <ul> + * <li> {@code lowerBound <= a < initial < b <= upperBound} </li> + * <li> {@code f(a) * f(b) <= 0} </li> + * </ul> + * If {@code f} is continuous on {@code [a,b]}, this means that {@code a} + * and {@code b} bracket a root of {@code f}. + * <p> + * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing + * values of k, where \( l_k = max(lower, initial - \delta_k) \), + * \( u_k = min(upper, initial + \delta_k) \), using recurrence + * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \). + * The algorithm stops when one of the following happens: <ul> + * <li> at least one positive and one negative value have been found -- success!</li> + * <li> both endpoints have reached their respective limits -- NoBracketingException </li> + * <li> {@code maximumIterations} iterations elapse -- NoBracketingException </li></ul> + * <p> + * If different signs are found at first iteration ({@code k=1}), then the returned + * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later + * iteration {@code k>1}, then the returned interval will be either + * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called + * with these parameters will therefore start with the smallest bracketing interval known + * at this step. + * </p> + * <p> + * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and + * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a + * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r} + * is larger than 1, the sequence has an asymptotically exponential rate. Note than the + * additive parameter {@code q} should never be set to zero, otherwise the interval would + * degenerate to the single initial point for all values of {@code k}. + * </p> + * <p> + * As a rule of thumb, when the location of the root is expected to be approximately known + * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the + * order of magnitude of the error margin. When the location of the root is really a wild guess, + * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval + * length at each iteration) and {@code q} should be set according to half the initial + * search interval length. + * </p> + * <p> + * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use + * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute + * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then + * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will + * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root. + * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned + * bracketing interval. + * </p> + * @param function function to check + * @param initial Initial midpoint of interval being expanded to + * bracket a root. + * @param lowerBound Lower bound (a is never lower than this value). + * @param upperBound Upper bound (b never is greater than this + * value). + * @param q additive offset used to compute bounds sequence (must be strictly positive) + * @param r multiplicative factor used to compute bounds sequence + * @param maximumIterations Maximum number of iterations to perform + * @return a two element array holding the bracketing values. + * @exception NoBracketingException if function cannot be bracketed in the search interval + */ + public static double[] bracket(final UnivariateFunction function, final double initial, + final double lowerBound, final double upperBound, + final double q, final double r, final int maximumIterations) + throws NoBracketingException { + + if (function == null) { + throw new NullArgumentException(LocalizedFormats.FUNCTION); + } + if (q <= 0) { + throw new NotStrictlyPositiveException(q); + } + if (maximumIterations <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations); + } + verifySequence(lowerBound, initial, upperBound); + + // initialize the recurrence + double a = initial; + double b = initial; + double fa = Double.NaN; + double fb = Double.NaN; + double delta = 0; + + for (int numIterations = 0; + (numIterations < maximumIterations) && (a > lowerBound || b < upperBound); + ++numIterations) { + + final double previousA = a; + final double previousFa = fa; + final double previousB = b; + final double previousFb = fb; + + delta = r * delta + q; + a = FastMath.max(initial - delta, lowerBound); + b = FastMath.min(initial + delta, upperBound); + fa = function.value(a); + fb = function.value(b); + + if (numIterations == 0) { + // at first iteration, we don't have a previous interval + // we simply compare both sides of the initial interval + if (fa * fb <= 0) { + // the first interval already brackets a root + return new double[] { a, b }; + } + } else { + // we have a previous interval with constant sign and expand it, + // we expect sign changes to occur at boundaries + if (fa * previousFa <= 0) { + // sign change detected at near lower bound + return new double[] { a, previousA }; + } else if (fb * previousFb <= 0) { + // sign change detected at near upper bound + return new double[] { previousB, b }; + } + } + + } + + // no bracketing found + throw new NoBracketingException(a, b, fa, fb); + + } + + /** + * Compute the midpoint of two values. + * + * @param a first value. + * @param b second value. + * @return the midpoint. + */ + public static double midpoint(double a, double b) { + return (a + b) * 0.5; + } + + /** + * Check whether the interval bounds bracket a root. That is, if the + * values at the endpoints are not equal to zero, then the function takes + * opposite signs at the endpoints. + * + * @param function Function. + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @return {@code true} if the function values have opposite signs at the + * given points. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static boolean isBracketing(UnivariateFunction function, + final double lower, + final double upper) + throws NullArgumentException { + if (function == null) { + throw new NullArgumentException(LocalizedFormats.FUNCTION); + } + final double fLo = function.value(lower); + final double fHi = function.value(upper); + return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0); + } + + /** + * Check whether the arguments form a (strictly) increasing sequence. + * + * @param start First number. + * @param mid Second number. + * @param end Third number. + * @return {@code true} if the arguments form an increasing sequence. + */ + public static boolean isSequence(final double start, + final double mid, + final double end) { + return (start < mid) && (mid < end); + } + + /** + * Check that the endpoints specify an interval. + * + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public static void verifyInterval(final double lower, + final double upper) + throws NumberIsTooLargeException { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, + lower, upper, false); + } + } + + /** + * Check that {@code lower < initial < upper}. + * + * @param lower Lower endpoint. + * @param initial Initial value. + * @param upper Upper endpoint. + * @throws NumberIsTooLargeException if {@code lower >= initial} or + * {@code initial >= upper}. + */ + public static void verifySequence(final double lower, + final double initial, + final double upper) + throws NumberIsTooLargeException { + verifyInterval(lower, initial); + verifyInterval(initial, upper); + } + + /** + * Check that the endpoints specify an interval and the end points + * bracket a root. + * + * @param function Function. + * @param lower Lower endpoint. + * @param upper Upper endpoint. + * @throws NoBracketingException if the function has the same sign at the + * endpoints. + * @throws NullArgumentException if {@code function} is {@code null}. + */ + public static void verifyBracketing(UnivariateFunction function, + final double lower, + final double upper) + throws NullArgumentException, + NoBracketingException { + if (function == null) { + throw new NullArgumentException(LocalizedFormats.FUNCTION); + } + verifyInterval(lower, upper); + if (!isBracketing(function, lower, upper)) { + throw new NoBracketingException(lower, upper, + function.value(lower), + function.value(upper)); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/solvers/package-info.java b/src/main/java/org/apache/commons/math3/analysis/solvers/package-info.java new file mode 100644 index 0000000..eb15fbc --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/solvers/package-info.java @@ -0,0 +1,22 @@ +/* + * 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. + */ +/** + * + * Root finding algorithms, for univariate real functions. + * + */ +package org.apache.commons.math3.analysis.solvers; |