| Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
| BrentSolver |
|
| 6.666666666666667;6.667 |
| 1 | /* |
|
| 2 | * Copyright 2003-2005 The Apache Software Foundation. |
|
| 3 | * |
|
| 4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
|
| 5 | * you may not use this file except in compliance with the License. |
|
| 6 | * You may obtain a copy of the License at |
|
| 7 | * |
|
| 8 | * http://www.apache.org/licenses/LICENSE-2.0 |
|
| 9 | * |
|
| 10 | * Unless required by applicable law or agreed to in writing, software |
|
| 11 | * distributed under the License is distributed on an "AS IS" BASIS, |
|
| 12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
| 13 | * See the License for the specific language governing permissions and |
|
| 14 | * limitations under the License. |
|
| 15 | */ |
|
| 16 | package org.apache.commons.math.analysis; |
|
| 17 | ||
| 18 | ||
| 19 | import org.apache.commons.math.ConvergenceException; |
|
| 20 | import org.apache.commons.math.FunctionEvaluationException; |
|
| 21 | ||
| 22 | /** |
|
| 23 | * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html"> |
|
| 24 | * Brent algorithm</a> for finding zeros of real univariate functions. |
|
| 25 | * <p> |
|
| 26 | * The function should be continuous but not necessarily smooth. |
|
| 27 | * |
|
| 28 | * @version $Revision$ $Date: 2005-06-03 22:36:42 -0700 (Fri, 03 Jun 2005) $ |
|
| 29 | */ |
|
| 30 | public class BrentSolver extends UnivariateRealSolverImpl { |
|
| 31 | ||
| 32 | /** Serializable version identifier */ |
|
| 33 | static final long serialVersionUID = 3350616277306882875L; |
|
| 34 | ||
| 35 | /** |
|
| 36 | * Construct a solver for the given function. |
|
| 37 | * |
|
| 38 | * @param f function to solve. |
|
| 39 | */ |
|
| 40 | public BrentSolver(UnivariateRealFunction f) { |
|
| 41 | 204 | super(f, 100, 1E-6); |
| 42 | 202 | } |
| 43 | ||
| 44 | /** |
|
| 45 | * Find a zero in the given interval. |
|
| 46 | * <p> |
|
| 47 | * Throws <code>ConvergenceException</code> if the values of the function |
|
| 48 | * at the endpoints of the interval have the same sign. |
|
| 49 | * |
|
| 50 | * @param min the lower bound for the interval. |
|
| 51 | * @param max the upper bound for the interval. |
|
| 52 | * @param initial the start value to use (ignored). |
|
| 53 | * @return the value where the function is zero |
|
| 54 | * @throws ConvergenceException the maximum iteration count is exceeded |
|
| 55 | * @throws FunctionEvaluationException if an error occurs evaluating |
|
| 56 | * the function |
|
| 57 | * @throws IllegalArgumentException if initial is not between min and max |
|
| 58 | */ |
|
| 59 | public double solve(double min, double max, double initial) |
|
| 60 | throws ConvergenceException, FunctionEvaluationException { |
|
| 61 | ||
| 62 | 0 | return solve(min, max); |
| 63 | } |
|
| 64 | ||
| 65 | /** |
|
| 66 | * Find a zero in the given interval. |
|
| 67 | * <p> |
|
| 68 | * Requires that the values of the function at the endpoints have opposite |
|
| 69 | * signs. An <code>IllegalArgumentException</code> is thrown if this is not |
|
| 70 | * the case. |
|
| 71 | * |
|
| 72 | * @param min the lower bound for the interval. |
|
| 73 | * @param max the upper bound for the interval. |
|
| 74 | * @return the value where the function is zero |
|
| 75 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
| 76 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
| 77 | * function |
|
| 78 | * @throws IllegalArgumentException if min is not less than max or the |
|
| 79 | * signs of the values of the function at the endpoints are not opposites |
|
| 80 | */ |
|
| 81 | public double solve(double min, double max) throws ConvergenceException, |
|
| 82 | FunctionEvaluationException { |
|
| 83 | ||
| 84 | 224 | clearResult(); |
| 85 | 224 | verifyInterval(min, max); |
| 86 | ||
| 87 | // Index 0 is the old approximation for the root. |
|
| 88 | // Index 1 is the last calculated approximation for the root. |
|
| 89 | // Index 2 is a bracket for the root with respect to x1. |
|
| 90 | 222 | double x0 = min; |
| 91 | 222 | double x1 = max; |
| 92 | double y0; |
|
| 93 | double y1; |
|
| 94 | 222 | y0 = f.value(x0); |
| 95 | 222 | y1 = f.value(x1); |
| 96 | ||
| 97 | // Verify bracketing |
|
| 98 | 222 | if (y0 * y1 >= 0) { |
| 99 | 8 | throw new IllegalArgumentException |
| 100 | ("Function values at endpoints do not have different signs." + |
|
| 101 | " Endpoints: [" + min + "," + max + "]" + |
|
| 102 | " Values: [" + y0 + "," + y1 + "]"); |
|
| 103 | } |
|
| 104 | ||
| 105 | 214 | double x2 = x0; |
| 106 | 214 | double y2 = y0; |
| 107 | 214 | double delta = x1 - x0; |
| 108 | 214 | double oldDelta = delta; |
| 109 | ||
| 110 | 214 | int i = 0; |
| 111 | 1602 | while (i < maximalIterationCount) { |
| 112 | 1602 | if (Math.abs(y2) < Math.abs(y1)) { |
| 113 | 404 | x0 = x1; |
| 114 | 404 | x1 = x2; |
| 115 | 404 | x2 = x0; |
| 116 | 404 | y0 = y1; |
| 117 | 404 | y1 = y2; |
| 118 | 404 | y2 = y0; |
| 119 | } |
|
| 120 | 1602 | if (Math.abs(y1) <= functionValueAccuracy) { |
| 121 | // Avoid division by very small values. Assume |
|
| 122 | // the iteration has converged (the problem may |
|
| 123 | // still be ill conditioned) |
|
| 124 | 12 | setResult(x1, i); |
| 125 | 12 | return result; |
| 126 | } |
|
| 127 | 1590 | double dx = (x2 - x1); |
| 128 | 1590 | double tolerance = |
| 129 | Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy); |
|
| 130 | 1590 | if (Math.abs(dx) <= tolerance) { |
| 131 | 202 | setResult(x1, i); |
| 132 | 202 | return result; |
| 133 | } |
|
| 134 | 1388 | if ((Math.abs(oldDelta) < tolerance) || |
| 135 | (Math.abs(y0) <= Math.abs(y1))) { |
|
| 136 | // Force bisection. |
|
| 137 | 18 | delta = 0.5 * dx; |
| 138 | 18 | oldDelta = delta; |
| 139 | } else { |
|
| 140 | 1370 | double r3 = y1 / y0; |
| 141 | double p; |
|
| 142 | double p1; |
|
| 143 | 1370 | if (x0 == x2) { |
| 144 | // Linear interpolation. |
|
| 145 | 924 | p = dx * r3; |
| 146 | 924 | p1 = 1.0 - r3; |
| 147 | } else { |
|
| 148 | // Inverse quadratic interpolation. |
|
| 149 | 446 | double r1 = y0 / y2; |
| 150 | 446 | double r2 = y1 / y2; |
| 151 | 446 | p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0)); |
| 152 | 446 | p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0); |
| 153 | } |
|
| 154 | 1370 | if (p > 0.0) { |
| 155 | 670 | p1 = -p1; |
| 156 | } else { |
|
| 157 | 700 | p = -p; |
| 158 | } |
|
| 159 | 1370 | if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) || |
| 160 | p >= Math.abs(0.5 * oldDelta * p1)) { |
|
| 161 | // Inverse quadratic interpolation gives a value |
|
| 162 | // in the wrong direction, or progress is slow. |
|
| 163 | // Fall back to bisection. |
|
| 164 | 64 | delta = 0.5 * dx; |
| 165 | 64 | oldDelta = delta; |
| 166 | } else { |
|
| 167 | 1306 | oldDelta = delta; |
| 168 | 1306 | delta = p / p1; |
| 169 | } |
|
| 170 | } |
|
| 171 | // Save old X1, Y1 |
|
| 172 | 1388 | x0 = x1; |
| 173 | 1388 | y0 = y1; |
| 174 | // Compute new X1, Y1 |
|
| 175 | 1388 | if (Math.abs(delta) > tolerance) { |
| 176 | 1172 | x1 = x1 + delta; |
| 177 | 216 | } else if (dx > 0.0) { |
| 178 | 122 | x1 = x1 + 0.5 * tolerance; |
| 179 | 94 | } else if (dx <= 0.0) { |
| 180 | 94 | x1 = x1 - 0.5 * tolerance; |
| 181 | } |
|
| 182 | 1388 | y1 = f.value(x1); |
| 183 | 1388 | if ((y1 > 0) == (y2 > 0)) { |
| 184 | 920 | x2 = x0; |
| 185 | 920 | y2 = y0; |
| 186 | 920 | delta = x1 - x0; |
| 187 | 920 | oldDelta = delta; |
| 188 | } |
|
| 189 | 1388 | i++; |
| 190 | } |
|
| 191 | 0 | throw new ConvergenceException("Maximum number of iterations exceeded."); |
| 192 | } |
|
| 193 | } |