# Brent's method in Julia

numerical-analysis root-finding julia

Brent's method or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines the bisection method, the secant method and inverse quadratic interpolation. This method always converges as long as the values of the function are computable within a given region containing a root.

Given a function $f(x)$ and the bracket $[x_0, x_1]$ two new points, $x_2$ and $x_3$, are initialized with the $x_1$ value. $x_0$ and $x_1$ are swapped if $|f(x_0)| < |f(x_1)|$.

Then, in each iteration if the evaluation of the points $x_0$, $x_1$ and $x_2$ are different (according to a certain tolerance) the inverse quadratic interpolation is used to get the new guess $x$. Otherwise, the linear interpolation (secant method) is used to obtain the guess. After that, if any of the following conditions are satisfied $x$ will be redefined using the bisection method:

• If $x$ does not lie between $(3x_0+x_1)/4$ and $x_1$.
• If in the previous iteration the bisection method was used or it is the first iteration and $|x-x_1| \geq {1 \over 2} \,|x_1-x_2|$.
• If in the previous iteration the bisection method was not used and $|x-x_1| \geq {1 \over 2} \,|x_2-x_3|$.
• If in the previous iteration the bisection method was used or it is the first iteration and $|x_1-x_2| < |\delta|$, where $\delta$ is a tolerance factor.
• If in the previous iteration the bisection method was not used and $|x_2-x_3| < |\delta|$, where $\delta$ is a tolerance factor.

We define $\delta$ as $2 \epsilon x_1$, where $\epsilon$ is the machine epsilon.

$x_3$ and $x_2$ are redefined in each iteration with $x_2$ and $x_1$ value, respectively, and the new guess $x$ will be set as $x_1$ if $f(x_0)f(x)<0$ or as $x_2$ otherwise. Again, $x_0$ and $x_1$ are swapped if $|f(x_0)| < |f(x_1)|$. The algorithm converges when $f(x)$ or $|x_1-x_0|$ are small enough, both according to tolerance factors.

In the following implementation, the inverse quadratic interpolation is applied directly. In other cases, like the implementation in Numerical recipes, used for example in Boost, the Lagrange polynomial is reduced defining the variables $p$, $q$, $r$, $s$ and $t$ as explained in MathWorld and $x$ value is not overwritten with the bisection method, but modified. For simplicity of the code, here the inverse quadratic interpolation is applied directly as in the entry Inverse quadratic interpolation in Julia and the new guess is overwritten if needed.

In Julia:

function brent(f::Function, x0::Number, x1::Number, args::Tuple=();
xtol::AbstractFloat=1e-7, ytol=2eps(Float64),
maxiter::Integer=50)
EPS = eps(Float64)
y0 = f(x0,args...)
y1 = f(x1,args...)
if abs(y0) < abs(y1)
# Swap lower and upper bounds.
x0, x1 = x1, x0
y0, y1 = y1, y0
end
x2 = x0
y2 = y0
x3 = x2
bisection = true
for _ in 1:maxiter
# x-tolerance.
if abs(x1-x0) < xtol
return x1
end

# Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
# and linear interpolation (secant method) otherwise.
if abs(y0-y2) > ytol && abs(y1-y2) > ytol
x = x0*y1*y2/((y0-y1)*(y0-y2)) +
x1*y0*y2/((y1-y0)*(y1-y2)) +
x2*y0*y1/((y2-y0)*(y2-y1))
else
x = x1 - y1 * (x1-x0)/(y1-y0)
end

# Use bisection method if satisfies the conditions.
delta = abs(2EPS*abs(x1))
min1 = abs(x-x1)
min2 = abs(x1-x2)
min3 = abs(x2-x3)
if (x < (3x0+x1)/4 && x > x1) ||
(bisection && min1 >= min2/2) ||
(!bisection && min1 >= min3/2) ||
(bisection && min2 < delta) ||
(!bisection && min3 < delta)
x = (x0+x1)/2
bisection = true
else
bisection = false
end

y = f(x,args...)
# y-tolerance.
if abs(y) < ytol
return x
end
x3 = x2
x2 = x1
if sign(y0) != sign(y)
x1 = x
y1 = y
else
x0 = x
y0 = y
end
if abs(y0) < abs(y1)
# Swap lower and upper bounds.
x0, x1 = x1, x0
y0, y1 = y1, y0
end
end
error("Max iteration exceeded")
end

As an example, for the function $f(x)=x^4-2x^2+1/4 \quad (0 \leq x \leq 1)$, the solution is $\sqrt{1-\sqrt{3}/2}$:

julia> (f(x)=x^4-2x^2+1/4; x0=0; x1=1; brent(f,x0,x1))
0.3660254037844386
julia> sqrt(1-sqrt(3)/2)
0.3660254037844387