What is "e" variable in popular implementations of Brent root finding algorithm?
I am reading the standard (Numerical Recipes and GSL C versions are identical) implementation of Brent root finding algorithm, and canno开发者_运维问答t understand the meaning of variable "e". The usage suggests that "e" is supposed to be the previous distance between the brackets. But then, why is it set to "xm" (half the distance) when we use bisection?
I'm not familiar with the algorithm. However, I can compare the C source and the Wikipedia description of the algorithm. The algorithm seems straight forward-ish (if you're familiar with methods to find roots), but the C implementation looks like a direct port of the fortran, so it's rather hard to read.
My best guess is that e
is related to the loop conditional.
Wikipedia says (line 8 of the algorithm):
repeat until f(b or s) = 0 or |b − a| is small enough (convergence)
The C source says:
e = b - a
, then later if (fabs(e) <= tol ...
.
I'd hope that the purpose of the variables would be described clearly in the book, but apparently not :)
Ok, here you go. I found the original implementation (in algol 60) here. In addition to a nice description of the algorithm, it says (starting on page 50):
let
e
be the value ofp/q
at the step before the last one. If|e|
< δ or|p/q|
≥1/2|e|
then do a bisection, otherwise we do either a bisection or interpolation just as in Dekker's algorithm. Thus|e|
decreases by at least a factor of two on every second step, and when|e|
< δ a bisection must be done. (After a bisection we takee = m
for the next step.)
So the addition of e
is Brent's "main modification" of Dekker's algorithm.
E is the "epsilon" variable, which is basically a measure of how close is close enough. Your particular application may not require 20 digits of precision, so epsilon lets you balance how many iterations it requires ( i.e., how long it runs ) versus how accurate you need it.
With floating point numbers you may not be able to be exact, so epsilon should be some small non-zero number. The actual value depends on your application... it's basically the largest acceptable error.
During a bisection step, the interval is exactly halved. Thus, e, holding the current width of the interval, is halved as well.
精彩评论