Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
SecantSolver |
|
| 5.0;5 |
1 | /* |
|
2 | * Copyright 2003-2004 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 | import java.io.Serializable; |
|
19 | ||
20 | import org.apache.commons.math.ConvergenceException; |
|
21 | import org.apache.commons.math.FunctionEvaluationException; |
|
22 | ||
23 | ||
24 | /** |
|
25 | * Implements a modified version of the |
|
26 | * <a href="http://mathworld.wolfram.com/SecantMethod.html">secant method</a> |
|
27 | * for approximating a zero of a real univariate function. |
|
28 | * <p> |
|
29 | * The algorithm is modified to maintain bracketing of a root by successive |
|
30 | * approximations. Because of forced bracketing, convergence may be slower than |
|
31 | * the unrestricted secant algorithm. However, this implementation should in |
|
32 | * general outperform the |
|
33 | * <a href="http://mathworld.wolfram.com/MethodofFalsePosition.html"> |
|
34 | * regula falsi method.</a> |
|
35 | * <p> |
|
36 | * The function is assumed to be continuous but not necessarily smooth. |
|
37 | * |
|
38 | * @version $Revision$ $Date: 2005-06-03 22:36:42 -0700 (Fri, 03 Jun 2005) $ |
|
39 | */ |
|
40 | public class SecantSolver extends UnivariateRealSolverImpl implements Serializable { |
|
41 | ||
42 | /** Serializable version identifier */ |
|
43 | static final long serialVersionUID = 1984971194738974867L; |
|
44 | ||
45 | /** |
|
46 | * Construct a solver for the given function. |
|
47 | * @param f function to solve. |
|
48 | */ |
|
49 | public SecantSolver(UnivariateRealFunction f) { |
|
50 | 8 | super(f, 100, 1E-6); |
51 | 6 | } |
52 | ||
53 | /** |
|
54 | * Find a zero in the given interval. |
|
55 | * |
|
56 | * @param min the lower bound for the interval |
|
57 | * @param max the upper bound for the interval |
|
58 | * @param initial the start value to use (ignored) |
|
59 | * @return the value where the function is zero |
|
60 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
61 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
62 | * function |
|
63 | * @throws IllegalArgumentException if min is not less than max or the |
|
64 | * signs of the values of the function at the endpoints are not opposites |
|
65 | */ |
|
66 | public double solve(double min, double max, double initial) |
|
67 | throws ConvergenceException, FunctionEvaluationException { |
|
68 | ||
69 | 0 | return solve(min, max); |
70 | } |
|
71 | ||
72 | /** |
|
73 | * Find a zero in the given interval. |
|
74 | * @param min the lower bound for the interval. |
|
75 | * @param max the upper bound for the interval. |
|
76 | * @return the value where the function is zero |
|
77 | * @throws ConvergenceException if the maximum iteration count is exceeded |
|
78 | * @throws FunctionEvaluationException if an error occurs evaluating the |
|
79 | * function |
|
80 | * @throws IllegalArgumentException if min is not less than max or the |
|
81 | * signs of the values of the function at the endpoints are not opposites |
|
82 | */ |
|
83 | public double solve(double min, double max) throws ConvergenceException, |
|
84 | FunctionEvaluationException { |
|
85 | ||
86 | 26 | clearResult(); |
87 | 26 | verifyInterval(min, max); |
88 | ||
89 | // Index 0 is the old approximation for the root. |
|
90 | // Index 1 is the last calculated approximation for the root. |
|
91 | // Index 2 is a bracket for the root with respect to x0. |
|
92 | // OldDelta is the length of the bracketing interval of the last |
|
93 | // iteration. |
|
94 | 26 | double x0 = min; |
95 | 26 | double x1 = max; |
96 | 26 | double y0 = f.value(x0); |
97 | 26 | double y1 = f.value(x1); |
98 | ||
99 | // Verify bracketing |
|
100 | 26 | if (y0 * y1 >= 0) { |
101 | 0 | throw new IllegalArgumentException |
102 | ("Function values at endpoints do not have different signs." + |
|
103 | " Endpoints: [" + min + "," + max + "]" + |
|
104 | " Values: [" + y0 + "," + y1 + "]"); |
|
105 | } |
|
106 | ||
107 | 26 | double x2 = x0; |
108 | 26 | double y2 = y0; |
109 | 26 | double oldDelta = x2 - x1; |
110 | 26 | int i = 0; |
111 | 216 | while (i < maximalIterationCount) { |
112 | 216 | if (Math.abs(y2) < Math.abs(y1)) { |
113 | 42 | x0 = x1; |
114 | 42 | x1 = x2; |
115 | 42 | x2 = x0; |
116 | 42 | y0 = y1; |
117 | 42 | y1 = y2; |
118 | 42 | y2 = y0; |
119 | } |
|
120 | 216 | if (Math.abs(y1) <= functionValueAccuracy) { |
121 | 10 | setResult(x1, i); |
122 | 10 | return result; |
123 | } |
|
124 | 206 | if (Math.abs(oldDelta) < |
125 | Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy)) { |
|
126 | 16 | setResult(x1, i); |
127 | 16 | return result; |
128 | } |
|
129 | double delta; |
|
130 | 190 | if (Math.abs(y1) > Math.abs(y0)) { |
131 | // Function value increased in last iteration. Force bisection. |
|
132 | 6 | delta = 0.5 * oldDelta; |
133 | } else { |
|
134 | 184 | delta = (x0 - x1) / (1 - y0 / y1); |
135 | 184 | if (delta / oldDelta > 1) { |
136 | // New approximation falls outside bracket. |
|
137 | // Fall back to bisection. |
|
138 | 4 | delta = 0.5 * oldDelta; |
139 | } |
|
140 | } |
|
141 | 190 | x0 = x1; |
142 | 190 | y0 = y1; |
143 | 190 | x1 = x1 + delta; |
144 | 190 | y1 = f.value(x1); |
145 | 190 | if ((y1 > 0) == (y2 > 0)) { |
146 | // New bracket is (x0,x1). |
|
147 | 126 | x2 = x0; |
148 | 126 | y2 = y0; |
149 | } |
|
150 | 190 | oldDelta = x2 - x1; |
151 | 190 | i++; |
152 | } |
|
153 | 0 | throw new ConvergenceException("Maximal iteration number exceeded" + i); |
154 | } |
|
155 | ||
156 | } |