CodeCogs

Calculates the zeros of a function using Brent's method.

Interface

#include <codecogs/maths/rootfinding/brent.h>

using namespace Maths::Rootfinding;

double brent (double (*f)(double), double x1 = -1E+7, double x2 = 1E+7, double eps = 1E-10, int maxit = 1000)

Brent's method is a root finding algorithm which combines root bracketing, interval bisection and inverse quadratic interpolation. It is sometimes known as the van Wijngaarden-Deker-Brent method.

The main algorithm uses a Lagrange interpolating polynomial of degree 2. Brent (1973) claims that this method will always converge as long as the values of the function are computable within a given region containing a root. Given three points $x_1$ , $x_2$ and $x_3$ , Brent's method fits x as a quadratic function of y, then uses the interpolation formula

$x = \sum_{k = 1}^3 x_k \prod_{1 \leq i \leq 3 \quad i \neq k} \frac {y - f(x_i)} {f(x_k) - f(x_i)}$

(1)

Subsequent root estimates are obtained by setting $y = 0$ , giving

$x = x_2 + \frac {P} {Q}$

(2)

where

$P = S[T(R - T)(x_3 - x_2) - (1 - R)(x_2 - x_1)] \f] \f[ Q = (T - 1)(R - 1)(S - 1)$

(3)

with

$R = \frac {f(x_2)} {f(x_3)} \qquad S = \frac {f(x_2)} {f(x_1)} \qquad T = \frac {f(x_1)} {f(x_3)}$

(4)

This algorithm finds the roots of the user-defined function f starting with an initial interval [x1, x2] and iterating the sequence above until either the accuracy eps is achieved or the maximum number of iterations maxit is exceeded.

References

Jean-Pierre Moreau's Home Page, http://perso.wanadoo.fr/jean-pierre.moreau/
MathWorld, http://mathworld.wolfram.com/BrentsMethod.html

Example 1

#include <codecogs/maths/rootfinding/brent.h>

#include <iostream>
#include <iomanip>

// user-defined function
double f(double x) {
    return (x + 1) * (x + 1) * (x + 1);
}

int main()  
{
  double x = Maths::RootFinding::brent(f, -3, 0);

  std::cout << "The calculated zero is X = " << std::setprecision(12) << x << std::endl;
  std::cout << "The associated ordinate value is Y = " << f(x) << std::endl;
  return 0;
}

Output:

The calculated zero is X = -0.999999878599
The associated ordinate value is Y = 1.78923630993e-021

Parameters

the user-defined function

Default value = -1E+7

Default value = 1E+7

eps

Default value = 1E-10

maxit

Default value = 1000