commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/solvers/BaseSecantSolver.java - commons-math - Git at Google

 /*
  * 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.math4.legacy.analysis.solvers;

 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
 import org.apache.commons.math4.legacy.exception.ConvergenceException;
 import org.apache.commons.math4.legacy.exception.MathInternalError;
 import org.apache.commons.math4.core.jdkmath.JdkMath;

 /**
  * 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} */
     @Override
     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} */
     @Override
     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 (JdkMath.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 (JdkMath.abs(x1 - x0) < JdkMath.max(rtol * JdkMath.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;

     }
 }
	/*
	* 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.math4.legacy.analysis.solvers;

	import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
	import org.apache.commons.math4.legacy.exception.ConvergenceException;
	import org.apache.commons.math4.legacy.exception.MathInternalError;
	import org.apache.commons.math4.core.jdkmath.JdkMath;

	/**
	* 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} */
	@Override
	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} */
	@Override
	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 (JdkMath.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 (JdkMath.abs(x1 - x0) < JdkMath.max(rtol * JdkMath.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;

	}
	}