1
2
3
4
5
6
7
8
9
10
11
12
13
14
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: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $
39 */
40 public class SecantSolver extends UnivariateRealSolverImpl implements Serializable {
41
42 /** Serializable version identifier */
43 private 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 super(f, 100, 1E-6);
51 }
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 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 clearResult();
87 verifyInterval(min, max);
88
89
90
91
92
93
94 double x0 = min;
95 double x1 = max;
96 double y0 = f.value(x0);
97 double y1 = f.value(x1);
98
99
100 if (y0 * y1 >= 0) {
101 throw new IllegalArgumentException
102 ("Function values at endpoints do not have different signs." +
103 " Endpoints: [" + min + "," + max + "]" +
104 " Values: [" + y0 + "," + y1 + "]");
105 }
106
107 double x2 = x0;
108 double y2 = y0;
109 double oldDelta = x2 - x1;
110 int i = 0;
111 while (i < maximalIterationCount) {
112 if (Math.abs(y2) < Math.abs(y1)) {
113 x0 = x1;
114 x1 = x2;
115 x2 = x0;
116 y0 = y1;
117 y1 = y2;
118 y2 = y0;
119 }
120 if (Math.abs(y1) <= functionValueAccuracy) {
121 setResult(x1, i);
122 return result;
123 }
124 if (Math.abs(oldDelta) <
125 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy)) {
126 setResult(x1, i);
127 return result;
128 }
129 double delta;
130 if (Math.abs(y1) > Math.abs(y0)) {
131
132 delta = 0.5 * oldDelta;
133 } else {
134 delta = (x0 - x1) / (1 - y0 / y1);
135 if (delta / oldDelta > 1) {
136
137
138 delta = 0.5 * oldDelta;
139 }
140 }
141 x0 = x1;
142 y0 = y1;
143 x1 = x1 + delta;
144 y1 = f.value(x1);
145 if ((y1 > 0) == (y2 > 0)) {
146
147 x2 = x0;
148 y2 = y0;
149 }
150 oldDelta = x2 - x1;
151 i++;
152 }
153 throw new ConvergenceException("Maximal iteration number exceeded" + i);
154 }
155
156 }