← Back to branch summary

~ubuntu-branches/ubuntu/natty/commons-math/natty

« back to all changes in this revision

Viewing changes to src/java/org/apache/commons/math/analysis/MullerSolver.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/java/org/apache/commons/math/analysis/MullerSolver.java
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;
18
 
 
19
 
import org.apache.commons.math.FunctionEvaluationException;
20
 
import org.apache.commons.math.MaxIterationsExceededException;
21
 
import org.apache.commons.math.util.MathUtils;
22
 
 
23
 
/**
24
 
 * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
25
 
 * Muller's Method</a> for root finding of real univariate functions. For
26
 
 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
27
 
 * chapter 3.
28
 
 * <p>
29
 
 * Muller's method applies to both real and complex functions, but here we
30
 
 * restrict ourselves to real functions. Methods solve() and solve2() find
31
 
 * real zeros, using different ways to bypass complex arithmetics.</p>
32
 
 *
33
 
 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
34
 
 * @since 1.2
35
 
 */
36
 
public class MullerSolver extends UnivariateRealSolverImpl {
37
 
 
38
 
    /** serializable version identifier */
39
 
    private static final long serialVersionUID = 6552227503458976920L;
40
 
 
41
 
    /**
42
 
     * Construct a solver for the given function.
43
 
     * 
44
 
     * @param f function to solve
45
 
     */
46
 
    public MullerSolver(UnivariateRealFunction f) {
47
 
        super(f, 100, 1E-6);
48
 
    }
49
 
 
50
 
    /**
51
 
     * Find a real root in the given interval with initial value.
52
 
     * <p>
53
 
     * Requires bracketing condition.</p>
54
 
     * 
55
 
     * @param min the lower bound for the interval
56
 
     * @param max the upper bound for the interval
57
 
     * @param initial the start value to use
58
 
     * @return the point at which the function value is zero
59
 
     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
60
 
     * or the solver detects convergence problems otherwise
61
 
     * @throws FunctionEvaluationException if an error occurs evaluating the
62
 
     * function
63
 
     * @throws IllegalArgumentException if any parameters are invalid
64
 
     */
65
 
    public double solve(double min, double max, double initial) throws
66
 
        MaxIterationsExceededException, FunctionEvaluationException {
67
 
 
68
 
        // check for zeros before verifying bracketing
69
 
        if (f.value(min) == 0.0) { return min; }
70
 
        if (f.value(max) == 0.0) { return max; }
71
 
        if (f.value(initial) == 0.0) { return initial; }
72
 
 
73
 
        verifyBracketing(min, max, f);
74
 
        verifySequence(min, initial, max);
75
 
        if (isBracketing(min, initial, f)) {
76
 
            return solve(min, initial);
77
 
        } else {
78
 
            return solve(initial, max);
79
 
        }
80
 
    }
81
 
 
82
 
    /**
83
 
     * Find a real root in the given interval.
84
 
     * <p>
85
 
     * Original Muller's method would have function evaluation at complex point.
86
 
     * Since our f(x) is real, we have to find ways to avoid that. Bracketing
87
 
     * condition is one way to go: by requiring bracketing in every iteration,
88
 
     * the newly computed approximation is guaranteed to be real.</p>
89
 
     * <p>
90
 
     * Normally Muller's method converges quadratically in the vicinity of a
91
 
     * zero, however it may be very slow in regions far away from zeros. For
92
 
     * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
93
 
     * bisection as a safety backup if it performs very poorly.</p>
94
 
     * <p>
95
 
     * The formulas here use divided differences directly.</p>
96
 
     * 
97
 
     * @param min the lower bound for the interval
98
 
     * @param max the upper bound for the interval
99
 
     * @return the point at which the function value is zero
100
 
     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
101
 
     * or the solver detects convergence problems otherwise
102
 
     * @throws FunctionEvaluationException if an error occurs evaluating the
103
 
     * function 
104
 
     * @throws IllegalArgumentException if any parameters are invalid
105
 
     */
106
 
    public double solve(double min, double max) throws MaxIterationsExceededException, 
107
 
        FunctionEvaluationException {
108
 
 
109
 
        // [x0, x2] is the bracketing interval in each iteration
110
 
        // x1 is the last approximation and an interpolation point in (x0, x2)
111
 
        // x is the new root approximation and new x1 for next round
112
 
        // d01, d12, d012 are divided differences
113
 
        double x0, x1, x2, x, oldx, y0, y1, y2, y;
114
 
        double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
115
 
 
116
 
        x0 = min; y0 = f.value(x0);
117
 
        x2 = max; y2 = f.value(x2);
118
 
        x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
119
 
 
120
 
        // check for zeros before verifying bracketing
121
 
        if (y0 == 0.0) { return min; }
122
 
        if (y2 == 0.0) { return max; }
123
 
        verifyBracketing(min, max, f);
124
 
 
125
 
        int i = 1;
126
 
        oldx = Double.POSITIVE_INFINITY;
127
 
        while (i <= maximalIterationCount) {
128
 
            // Muller's method employs quadratic interpolation through
129
 
            // x0, x1, x2 and x is the zero of the interpolating parabola.
130
 
            // Due to bracketing condition, this parabola must have two
131
 
            // real roots and we choose one in [x0, x2] to be x.
132
 
            d01 = (y1 - y0) / (x1 - x0);
133
 
            d12 = (y2 - y1) / (x2 - x1);
134
 
            d012 = (d12 - d01) / (x2 - x0);
135
 
            c1 = d01 + (x1 - x0) * d012;
136
 
            delta = c1 * c1 - 4 * y1 * d012;
137
 
            xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
138
 
            xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
139
 
            // xplus and xminus are two roots of parabola and at least
140
 
            // one of them should lie in (x0, x2)
141
 
            x = isSequence(x0, xplus, x2) ? xplus : xminus;
142
 
            y = f.value(x);
143
 
 
144
 
            // check for convergence
145
 
            tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
146
 
            if (Math.abs(x - oldx) <= tolerance) {
147
 
                setResult(x, i);
148
 
                return result;
149
 
            }
150
 
            if (Math.abs(y) <= functionValueAccuracy) {
151
 
                setResult(x, i);
152
 
                return result;
153
 
            }
154
 
 
155
 
            // Bisect if convergence is too slow. Bisection would waste
156
 
            // our calculation of x, hopefully it won't happen often.
157
 
            // the real number equality test x == x1 is intentional and
158
 
            // completes the proximity tests above it
159
 
            boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
160
 
                             (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
161
 
                             (x == x1);
162
 
            // prepare the new bracketing interval for next iteration
163
 
            if (!bisect) {
164
 
                x0 = x < x1 ? x0 : x1;
165
 
                y0 = x < x1 ? y0 : y1;
166
 
                x2 = x > x1 ? x2 : x1;
167
 
                y2 = x > x1 ? y2 : y1;
168
 
                x1 = x; y1 = y;
169
 
                oldx = x;
170
 
            } else {
171
 
                double xm = 0.5 * (x0 + x2);
172
 
                double ym = f.value(xm);
173
 
                if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
174
 
                    x2 = xm; y2 = ym;
175
 
                } else {
176
 
                    x0 = xm; y0 = ym;
177
 
                }
178
 
                x1 = 0.5 * (x0 + x2);
179
 
                y1 = f.value(x1);
180
 
                oldx = Double.POSITIVE_INFINITY;
181
 
            }
182
 
            i++;
183
 
        }
184
 
        throw new MaxIterationsExceededException(maximalIterationCount);
185
 
    }
186
 
 
187
 
    /**
188
 
     * Find a real root in the given interval.
189
 
     * <p>
190
 
     * solve2() differs from solve() in the way it avoids complex operations.
191
 
     * Except for the initial [min, max], solve2() does not require bracketing
192
 
     * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
193
 
     * number arises in the computation, we simply use its modulus as real
194
 
     * approximation.</p>
195
 
     * <p>
196
 
     * Because the interval may not be bracketing, bisection alternative is
197
 
     * not applicable here. However in practice our treatment usually works
198
 
     * well, especially near real zeros where the imaginary part of complex
199
 
     * approximation is often negligible.</p>
200
 
     * <p>
201
 
     * The formulas here do not use divided differences directly.</p>
202
 
     * 
203
 
     * @param min the lower bound for the interval
204
 
     * @param max the upper bound for the interval
205
 
     * @return the point at which the function value is zero
206
 
     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
207
 
     * or the solver detects convergence problems otherwise
208
 
     * @throws FunctionEvaluationException if an error occurs evaluating the
209
 
     * function 
210
 
     * @throws IllegalArgumentException if any parameters are invalid
211
 
     */
212
 
    public double solve2(double min, double max) throws MaxIterationsExceededException, 
213
 
        FunctionEvaluationException {
214
 
 
215
 
        // x2 is the last root approximation
216
 
        // x is the new approximation and new x2 for next round
217
 
        // x0 < x1 < x2 does not hold here
218
 
        double x0, x1, x2, x, oldx, y0, y1, y2, y;
219
 
        double q, A, B, C, delta, denominator, tolerance;
220
 
 
221
 
        x0 = min; y0 = f.value(x0);
222
 
        x1 = max; y1 = f.value(x1);
223
 
        x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
224
 
 
225
 
        // check for zeros before verifying bracketing
226
 
        if (y0 == 0.0) { return min; }
227
 
        if (y1 == 0.0) { return max; }
228
 
        verifyBracketing(min, max, f);
229
 
 
230
 
        int i = 1;
231
 
        oldx = Double.POSITIVE_INFINITY;
232
 
        while (i <= maximalIterationCount) {
233
 
            // quadratic interpolation through x0, x1, x2
234
 
            q = (x2 - x1) / (x1 - x0);
235
 
            A = q * (y2 - (1 + q) * y1 + q * y0);
236
 
            B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
237
 
            C = (1 + q) * y2;
238
 
            delta = B * B - 4 * A * C;
239
 
            if (delta >= 0.0) {
240
 
                // choose a denominator larger in magnitude
241
 
                double dplus = B + Math.sqrt(delta);
242
 
                double dminus = B - Math.sqrt(delta);
243
 
                denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
244
 
            } else {
245
 
                // take the modulus of (B +/- Math.sqrt(delta))
246
 
                denominator = Math.sqrt(B * B - delta);
247
 
            }
248
 
            if (denominator != 0) {
249
 
                x = x2 - 2.0 * C * (x2 - x1) / denominator;
250
 
                // perturb x if it exactly coincides with x1 or x2
251
 
                // the equality tests here are intentional
252
 
                while (x == x1 || x == x2) {
253
 
                    x += absoluteAccuracy;
254
 
                }
255
 
            } else {
256
 
                // extremely rare case, get a random number to skip it
257
 
                x = min + Math.random() * (max - min);
258
 
                oldx = Double.POSITIVE_INFINITY;
259
 
            }
260
 
            y = f.value(x);
261
 
 
262
 
            // check for convergence
263
 
            tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
264
 
            if (Math.abs(x - oldx) <= tolerance) {
265
 
                setResult(x, i);
266
 
                return result;
267
 
            }
268
 
            if (Math.abs(y) <= functionValueAccuracy) {
269
 
                setResult(x, i);
270
 
                return result;
271
 
            }
272
 
 
273
 
            // prepare the next iteration
274
 
            x0 = x1; y0 = y1;
275
 
            x1 = x2; y1 = y2;
276
 
            x2 = x; y2 = y;
277
 
            oldx = x;
278
 
            i++;
279
 
        }
280
 
        throw new MaxIterationsExceededException(maximalIterationCount);
281
 
    }
282
 
}