Rev 3408 | Blame | Compare with Previous | Last modification | View Log | RSS feed
/*
* Copyright (C) 2009 Christian Gawron
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
*
* Author: Christian Gawron
* Create date: 03-Jul-2009
*/
package uk.me.parabola.mkgmap.reader.dem;
/**
* Find zero of a function using Brent's method.
*/
public class Brent {
public interface Function {
public double eval(double x);
}
private static final double epsilon = 3.0e-10;
public static double zero(Function f, double x1, double x2) {
return zero(f, x1, x2, 1e-8, 100);
}
public static double zero(Function f, double x1, double x2, double tol, int maxit) {
double a=x1, b=x2, c=x2;
double fa=f.eval(a), fb=f.eval(b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
throw new ArithmeticException("Root must be bracketed");
double fc = fb;
double d = 0;
double e = 0;
for (int iterations=0; iterations < maxit; iterations++) {
if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
c=a;
fc=fa;
e=d=b-a;
}
if (Math.abs(fc) < Math.abs(fb)) {
a=b;
b=c;
c=a;
fa=fb;
fb=fc;
fc=fa;
}
double tolerance = 2.0 * epsilon * Math.abs(b) + 0.5 * tol;
double xm = 0.5 * (c - b);
if (Math.abs(xm) <= tolerance || fb == 0.0) return b;
if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) {
double s = fb / fa;
double p;
double q;
if (a == c) {
p=2.0*xm*s;
q=1.0-s;
} else {
q=fa/fc;
double r = fb / fc;
p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p > 0.0) q = -q;
p=Math.abs(p);
double min1 = 3.0 * xm * q - Math.abs(tolerance * q);
double min2 = Math.abs(e * q);
if (2.0*p < (min1 < min2 ? min1 : min2)) {
e=d;
d=p/q;
} else {
d=xm;
e=d;
}
} else {
d=xm;
e=d;
}
a=b;
fa=fb;
if (Math.abs(d) > tolerance)
b += d;
else
b += xm >= 0 ? tolerance : -tolerance;
fb=f.eval(b);
}
throw new ArithmeticException("Maximum number of iterations exceeded");
}
}