1 /*
2 * jcurl java curling software framework http://www.jcurl.org Copyright (C)
3 * 2005-2009 M. Rohrmoser
4 *
5 * This program is free software; you can redistribute it and/or modify it under
6 * the terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 2 of the License, or (at your option) any later
8 * version.
9 *
10 * This program is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License along with
16 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple
17 * Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19
20 package org.jcurl.math;
21
22 import org.apache.commons.logging.Log;
23 import org.jcurl.core.log.JCLoggerFactory;
24
25 /**
26 * Bisection root-finding algorithm.
27 *
28 * @author <a href="mailto:m@jcurl.org">M. Rohrmoser </a>
29 * @version $Id$
30 */
31 public class BisectionSolver {
32 private static final Log log = JCLoggerFactory
33 .getLogger(BisectionSolver.class);
34
35 /**
36 * Bisection root-finding algorithm. Ported from
37 * http://en.wikipedia.org/wiki/Bisection_method
38 *
39 * @param f
40 * the function to test
41 * @param y
42 * y value to find x for. 0 finds the root.
43 * @param left
44 * lower end of the interval to scan
45 * @param right
46 * upper end of the interval to scan
47 * @param epsilon
48 * iteration stops when
49 * <code>Math.abs(right - left) <= epsilon</code>
50 */
51 public static double findRoot(final R1R1Function f, final double y,
52 double left, double right, final double epsilon) {
53 double lefty = f.at(left, 0);
54 if (Double.isNaN(lefty))
55 return Double.NaN;
56 // check if left-y an right-y have opposite signs. Use signum to avoid
57 // overflows.
58 if (Math.signum(lefty) * f.at(right, 0) > 0)
59 return Double.NaN;
60 double midy = Double.NaN;
61 while (Math.abs(right - left) > epsilon) {
62 // Calculate midpoint of domain
63 final double midpoint = (right + left) / 2;
64 midy = f.at(midpoint, 0);
65 if (log.isDebugEnabled())
66 log.debug("f(" + left + ")=" + lefty + " f(" + midpoint + ")="
67 + midy + " f(" + right + ")=unused");
68 if (Double.isNaN(midy))
69 return Double.NaN;
70 // examine midy
71 if ((lefty - y) * (midy - y) > 0) {
72 // Throw away left half
73 left = midpoint;
74 lefty = midy;
75 } else
76 // Throw away right half
77 right = midpoint;
78 }
79 return (right + left) / 2;
80 }
81 }