View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.analysis.solvers;
18  
19  import org.apache.commons.math.ConvergenceException;
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MaxIterationsExceededException;
22  import org.apache.commons.math.analysis.UnivariateRealFunction;
23  import org.apache.commons.math.util.MathUtils;
24  
25  /**
26   * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
27   * Ridders' Method</a> for root finding of real univariate functions. For
28   * reference, see C. Ridders, <i>A new algorithm for computing a single root
29   * of a real continuous function </i>, IEEE Transactions on Circuits and
30   * Systems, 26 (1979), 979 - 980.
31   * <p>
32   * The function should be continuous but not necessarily smooth.</p>
33   *  
34   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
35   * @since 1.2
36   */
37  public class RiddersSolver extends UnivariateRealSolverImpl {
38  
39      /**
40       * Construct a solver for the given function.
41       * 
42       * @param f function to solve
43       * @deprecated as of 2.0 the function to solve is passed as an argument
44       * to the {@link #solve(UnivariateRealFunction, double, double)} or
45       * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
46       * method.
47       */
48      @Deprecated
49      public RiddersSolver(UnivariateRealFunction f) {
50          super(f, 100, 1E-6);
51      }
52  
53      /**
54       * Construct a solver.
55       */
56      public RiddersSolver() {
57          super(100, 1E-6);
58      }
59  
60      /** {@inheritDoc} */
61      @Deprecated
62      public double solve(final double min, final double max)
63          throws ConvergenceException, FunctionEvaluationException {
64          return solve(f, min, max);
65      }
66  
67      /** {@inheritDoc} */
68      @Deprecated
69      public double solve(final double min, final double max, final double initial)
70          throws ConvergenceException, FunctionEvaluationException {
71          return solve(f, min, max, initial);
72      }
73  
74      /**
75       * Find a root in the given interval with initial value.
76       * <p>
77       * Requires bracketing condition.</p>
78       * 
79       * @param f the function to solve
80       * @param min the lower bound for the interval
81       * @param max the upper bound for the interval
82       * @param initial the start value to use
83       * @return the point at which the function value is zero
84       * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
85       * @throws FunctionEvaluationException if an error occurs evaluating the
86       * function
87       * @throws IllegalArgumentException if any parameters are invalid
88       */
89      public double solve(final UnivariateRealFunction f,
90                          final double min, final double max, final double initial)
91          throws MaxIterationsExceededException, FunctionEvaluationException {
92  
93          // check for zeros before verifying bracketing
94          if (f.value(min) == 0.0) { return min; }
95          if (f.value(max) == 0.0) { return max; }
96          if (f.value(initial) == 0.0) { return initial; }
97  
98          verifyBracketing(min, max, f);
99          verifySequence(min, initial, max);
100         if (isBracketing(min, initial, f)) {
101             return solve(f, min, initial);
102         } else {
103             return solve(f, initial, max);
104         }
105     }
106 
107     /**
108      * Find a root in the given interval.
109      * <p>
110      * Requires bracketing condition.</p>
111      * 
112      * @param f the function to solve
113      * @param min the lower bound for the interval
114      * @param max the upper bound for the interval
115      * @return the point at which the function value is zero
116      * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
117      * @throws FunctionEvaluationException if an error occurs evaluating the
118      * function 
119      * @throws IllegalArgumentException if any parameters are invalid
120      */
121     public double solve(final UnivariateRealFunction f,
122                         final double min, final double max)
123         throws MaxIterationsExceededException, FunctionEvaluationException {
124 
125         // [x1, x2] is the bracketing interval in each iteration
126         // x3 is the midpoint of [x1, x2]
127         // x is the new root approximation and an endpoint of the new interval
128         double x1, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance;
129 
130         x1 = min; y1 = f.value(x1);
131         x2 = max; y2 = f.value(x2);
132 
133         // check for zeros before verifying bracketing
134         if (y1 == 0.0) { return min; }
135         if (y2 == 0.0) { return max; }
136         verifyBracketing(min, max, f);
137 
138         int i = 1;
139         oldx = Double.POSITIVE_INFINITY;
140         while (i <= maximalIterationCount) {
141             // calculate the new root approximation
142             x3 = 0.5 * (x1 + x2);
143             y3 = f.value(x3);
144             if (Math.abs(y3) <= functionValueAccuracy) {
145                 setResult(x3, i);
146                 return result;
147             }
148             delta = 1 - (y1 * y2) / (y3 * y3);  // delta > 1 due to bracketing
149             correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
150                          (x3 - x1) / Math.sqrt(delta);
151             x = x3 - correction;                // correction != 0
152             y = f.value(x);
153 
154             // check for convergence
155             tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
156             if (Math.abs(x - oldx) <= tolerance) {
157                 setResult(x, i);
158                 return result;
159             }
160             if (Math.abs(y) <= functionValueAccuracy) {
161                 setResult(x, i);
162                 return result;
163             }
164 
165             // prepare the new interval for next iteration
166             // Ridders' method guarantees x1 < x < x2
167             if (correction > 0.0) {             // x1 < x < x3
168                 if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
169                     x2 = x; y2 = y;
170                 } else {
171                     x1 = x; x2 = x3;
172                     y1 = y; y2 = y3;
173                 }
174             } else {                            // x3 < x < x2
175                 if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
176                     x1 = x; y1 = y;
177                 } else {
178                     x1 = x3; x2 = x;
179                     y1 = y3; y2 = y;
180                 }
181             }
182             oldx = x;
183             i++;
184         }
185         throw new MaxIterationsExceededException(maximalIterationCount);
186     }
187 }