Ctrl+K

Root finding

9. Root finding#

9.1. One-dimensional root finding#

Often times in computational settings, we wish to find the roots of a system. Some simple examples are find the roots of a denominator, to see where poles arise, and avoid them during your numerical integation scheme. The two most common root finding algorithms are Newton’s Method and bisection method. The former has a more numerically stable counter part, calle the secant method, while the secant method and bisection method have been combined to form the Brent method. More optimized formulas exist for root finding, with very specific use cases. However, for this chaper, we will focus on the secant, bisection and Brent methods.

9.1.1. Secant method#

Given a function f(x), and an initial guesses x0 and x1 for the location of the root x∗, the secant method iterates the formula

xn=xn−1−f(xn−1)xn−1−xn−2f(xn−1)−f(xn−1).

Then, $lim_{n \to \infty} x_n = x_\ast.

In code, this can be implemented very simply

// CODE HAS NOT BEEN CHECKED
template<typename T, typename Func, typename...Args>
std::optional<T>
secant_method(Func&& func, T x_0, T x_1, T tol, int max_iter, Args&&... args)
{
    T x_n;
    do 
    {
        T f_n1 = func(x_1, std::forward<Args>(args)...);
        f f_n0 = func(x_0, std::forward<Args>(args)...);
        x_n = x_1 - f_n1 * (x_n1 - x_n0) / (f_n1 - f_n0);

        if (std::fabs(func(x_n, std::forward<Args>(args)...) < tol) return x_n;
        else {
            x_0 = x_1;
            x_1 = x_n;
        }
        max_iter--;
    } while (max_iter > 0);

    return std::nullopt;
}

As it is posible for the method to not converge within the allotted ierations, we provide a safety mechanism for the user the test during run time, through returning an optional. That’s it!

9.1.2. Bisection method#

As the name suggests, this method finds the root by systematically splitting an interval into two, and continuing the search in one of those intervals. Let a and b be the end points of the interval, such that a<b, and f(x) be the function, whose roots we are locating. If f(a)f(b)>0, we can conclude tha there is no sign change for f(x), and consequently (though not strictly true), there is no root. That being said, you should always choose the initial interval, such that f(a)f(b)<0.

The bisection method uses the same logic as just described: let c0 be the midpoint between a and b, if f(a)f(c0)<0, the root is in the first half of the interval, else it most be in the other half. Let’s assumme the product is negative, then we let c1 be the midpoint between a and c0 and check if f(a)f(c1)<0, and so on. Written in code, the bisection method can be implemented as

// CODE HAS NOT BEEN CHECKED
template<typename T, typename Func, typename...Args>
std::optional<T>
bisection_method(Func&& func, T a, T b, T tol, int& max_iter, Args&&... args)
{
    // Note that a empty optional is now ambiguous as to the error.
    // However, by passint `max_iter` as a reference, we can check if it's zero
    if (f(a) * f(b) > 0) return std::nullopt; 
    if (a > b) std::swap(a, b);

    T c{ (b - a) / 2.0 };
    do {
        if (std::fabs(c) < tol) return c;
        
        if (f(a) * f(c) < 0) {
            c = b;
        } else {
            c = a;
        }
    } while(max_iter > 0);

    return std::nullopt;
}

9.1.3. Brent method#

Contents