org.apache.commons.math.analysis.solvers.RiddersSolver Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-math Show documentation
Show all versions of commons-math Show documentation
The Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.
/*
* 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.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the
* Ridders' Method for root finding of real univariate functions. For
* reference, see C. Ridders, A new algorithm for computing a single root
* of a real continuous function , IEEE Transactions on Circuits and
* Systems, 26 (1979), 979 - 980.
*
* The function should be continuous but not necessarily smooth.
*
* @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
* @since 1.2
*/
public class RiddersSolver extends UnivariateRealSolverImpl {
/**
* Construct a solver for the given function.
*
* @param f function to solve
* @deprecated as of 2.0 the function to solve is passed as an argument
* to the {@link #solve(UnivariateRealFunction, double, double)} or
* {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
* method.
*/
@Deprecated
public RiddersSolver(UnivariateRealFunction f) {
super(f, 100, 1E-6);
}
/**
* Construct a solver.
*/
public RiddersSolver() {
super(100, 1E-6);
}
/** {@inheritDoc} */
@Deprecated
public double solve(final double min, final double max)
throws ConvergenceException, FunctionEvaluationException {
return solve(f, min, max);
}
/** {@inheritDoc} */
@Deprecated
public double solve(final double min, final double max, final double initial)
throws ConvergenceException, FunctionEvaluationException {
return solve(f, min, max, initial);
}
/**
* Find a root in the given interval with initial value.
*
* Requires bracketing condition.
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws MaxIterationsExceededException, FunctionEvaluationException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(f, min, initial);
} else {
return solve(f, initial, max);
}
}
/**
* Find a root in the given interval.
*
* Requires bracketing condition.
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, FunctionEvaluationException {
// [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, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance;
x1 = min; y1 = f.value(x1);
x2 = max; y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y1 == 0.0) { return min; }
if (y2 == 0.0) { return max; }
verifyBracketing(min, max, f);
int i = 1;
oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
// calculate the new root approximation
x3 = 0.5 * (x1 + x2);
y3 = f.value(x3);
if (Math.abs(y3) <= functionValueAccuracy) {
setResult(x3, i);
return result;
}
delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
(x3 - x1) / Math.sqrt(delta);
x = x3 - correction; // correction != 0
y = f.value(x);
// check for convergence
tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
if (Math.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (Math.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the new interval for next iteration
// Ridders' method guarantees x1 < x < x2
if (correction > 0.0) { // x1 < x < x3
if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
x2 = x; y2 = y;
} else {
x1 = x; x2 = x3;
y1 = y; y2 = y3;
}
} else { // x3 < x < x2
if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
x1 = x; y1 = y;
} else {
x1 = x3; x2 = x;
y1 = y3; y2 = y;
}
}
oldx = x;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}