Ridder's Algorithm
Introduction
Ridder’s algorithm is a rootfinding algorithm that helps to find the root of a given function. This algorithm is an improvement over the bisection method and provides a faster convergence rate. The algorithm is named after Derek Ridder, who introduced it in 1978. It is widely used in scientific and engineering applications.
Implementation
The algorithm is implemented in OCaml as follows:
let ridder ?(max_iter = 1000) ?(xtol = 1e6) f a b =
let fa = f a in
let fb = f b in
let error () =
let s = Printf.sprintf "f(a) *. f(b) = %g *. %g should be negative." fa fb in
Owl_exception.INVALID_ARGUMENT s
in
Owl_exception.verify (fa *. fb < 0.) error;
if fa = 0.
then a
else if fb = 0.
then b
else (
let xa = ref a in
let xb = ref b in
let fa = ref fa in
let fb = ref fb in
let x = ref infinity in
try
for _i = 1 to max_iter do
let dm = 0.5 *. (!xb . !xa) in
let xm = !xa +. dm in
let fm = f xm in
let s = sqrt ((fm *. fm) . (!fa *. !fb)) in
assert (s <> 0.);
let sgn = if !fa < !fb then 1. else 1. in
x := xm +. (sgn *. dm *. fm /. s);
let fn = f !x in
if Owl_base_maths.same_sign fn fm = false
then (
xa := !x;
xb := xm;
fa := fn;
fb := fm)
else if Owl_base_maths.same_sign fn !fa = false
then (
xb := !x;
fb := fn)
else (
xa := !x;
fa := fn);
assert (abs_float (!xb . !xa) >= xtol && fn != 0.)
done;
!x
with
 _ > !x)
The function takes three arguments: the function f
whose root needs to be found, and the interval [a, b]
within which to search for the root. The function returns the estimated value of the root.
Stepbystep Explanation

The function
ridder
takes three input parameters: the functionf
whose root needs to be found, and the interval[a, b]
within which to search for the root. It also has two optional parameters:max_iter
(maximum number of iterations) andxtol
(tolerance level). The function returns the estimated value of the root. 
The function computes the value of
f
at the two endpointsa
andb
, and stores them in the variablesfa
andfb
, respectively. 
The function checks if the product of
fa
andfb
is negative. If it is not, it raises an error since the algorithm only works if the function changes sign over the interval[a, b]
. 
The function checks if either
fa
orfb
is zero. If either of them is zero, the function returns the corresponding value ofa
orb
since the root is already found. 
The function initializes the variables
xa
,xb
,fa
, andfb
with the values ofa
,b
,fa
, andfb
, respectively. It also initializes the variablex
with a value of infinity. 
The function enters a loop that runs for a maximum of
max_iter
iterations. This is to prevent the function from running indefinitely in case the algorithm fails to converge. 
Within the loop, the function computes the midpoint
xm
of the interval[xa, xb]
. 
The function computes the value of
f
atxm
and stores it in the variablefm
. 
The function computes the value of
s
, which is the square root of(fm^2  fa*fb)
. This value is used to determine the new estimate for the root. 
The function checks if
s
is zero. If it is, the algorithm cannot continue since it would require division by zero. In this case, the function raises an error. 
The function computes the sign of
fa
andfb
and stores it in the variablesgn
. Iffa < fb
,sgn
is set to1
, otherwise it is set to1
. 
The function computes the new estimate for the root using the formula
x = xm + (sgn * dm * fm / s)
, wheredm
is half the distance betweenxa
andxb
. 
The function computes the value of
f
at the new estimatex
and stores it in the variablefn
. 
The function checks if
fn
andfm
have opposite signs. If they do, the root is betweenxm
andx
, so the interval is updated to[xa, x]
andfa
andfb
are updated accordingly. If they don’t, the root is betweenxa
andx
, so the interval is updated to[x, xb]
andfa
andfb
are also updated accordingly. 
The function checks if
fn
andfa
have opposite signs. If they do, the root is betweenxa
andx
, so the interval is updated to[xa, x]
andfa
is updated tofn
. If they don’t, the root is betweenx
andxb
, so the interval is updated to[x, xb]
andfb
is updated tofn
. 
The function checks if the absolute difference between
xa
andxb
is less thanxtol
. If it is, the algorithm has converged and the function returns the estimated value of the root. 
If the loop has completed
max_iter
iterations without converging, the function raises an error. 
If an error occurs during the loop, the function returns the current estimate of the root.
Overall, the algorithm works by iteratively refining the interval [xa, xb]
until the root is found within a certain tolerance level. At each iteration, the algorithm computes the midpoint xm
of the interval and the value of f
at xm
. It then uses these values to compute a new estimate for the root using the Ridder’s formula. The algorithm then updates the interval and the values of fa
and fb
based on whether the new estimate is closer to xa
or xb
. The algorithm repeats this process until the root is found or the maximum number of iterations is reached.
Complexity Analysis
The time complexity of the Ridder’s algorithm is O(log2(1/ε))
, where ε
is the tolerance level. This means that the algorithm converges exponentially fast as the tolerance level is decreased. The algorithm also has a space complexity of O(1)
since it only uses a fixed number of variables regardless of the size of the input.
However, it is important to note that the algorithm may not converge if the function is not wellbehaved or if the interval [a, b]
is not chosen properly. In such cases, the algorithm may require more iterations or fail to converge altogether. Therefore, it is important to carefully choose the interval and verify that the function changes sign over the interval before using the algorithm.