One-Dimensional Root-Finding

This section delves into a range of classes and functions dedicated to the finding roots in one-dimensional functions. The featured algorithms are divided into two main types: bracketing algorithms, which operate without the need for function derivatives, and polishing algorithms, which require the computation of the function’s derivative. In addition, the library offers search algorithms that can be used to find a bracket where a root exists when the initial guess is not near the actual root.

All functions are available through the Roots.hpp header file.

Overview

Popular one-dimensional root-finding methods encompass the bisection method, Newton’s method, Ridder’s method, and the secant method. Each of these methods brings its unique advantages and potential drawbacks, and their applicability hinges on the nature of the problem at hand.

Bracketing methods such as the bisection method and Ridder’s method are designed to locate a root within a specified interval. These methods are particularly useful when the function’s derivative is unavailable or when the function is not differentiable. The bisection method, in particular, is known for its simplicity and reliability, as it guarantees convergence to a root as long as the function is continuous and changes sign within the interval. Ridder’s method, on the other hand, can converge with fewer iterations, but each iteration is more computationally expensive.

The following bracketing methods are available in the library:

  • Bisection Method: Known for its reliability and surefire convergence, the bisection method is a good choice when you need a method that is guaranteed to find a root. However, it can sometimes converge at a slower pace compared to other methods. (Wikipedia: Bisection Method)

  • Ridder’s Method: Outpacing the bisection method, Ridder’s method can converge more rapidly. However, each iteration is more computationally expensive, and the function must have a continuous second derivative. (Wikipedia: Ridder’s Method)

  • Regula Falsi Method: Also known as the false position method, the regula falsi method is a bracketing algorithm that uses linear interpolation to find the root of a function. It is similar to the bisection method, but it can converge faster in some cases. (Wikipedia: Regula Falsi Method)

Polishing methods, such as Newton’s method and Steffensen’s method, are designed to refine an initial guess to find a root. These methods are generally faster than bracketing methods, but they require the first derivative of the function to be available. Newton’s method, in particular, is known for its speed and efficiency, but it may fail to converge or converge to a local minimum instead of a root.

The following polishing methods are available in the library:

  • Newton’s Method: A popular root-finding algorithm, Newton’s method is known for its speed and efficiency. However, it requires the function to be differentiable and the derivative to be non-zero at the estimate. (Wikipedia: Newton’s Method)

  • Secant Method: The secant method is a root-finding algorithm that does not require the derivative of the function. It is similar to the Newton-Raphson method, but it uses a finite difference approximation to the derivative. (Wikipedia: Secant Method)

  • Steffensen’s Method: Steffensen’s method is an iterative root-finding algorithm that improves upon the simple fixed-point iteration by incorporating a form of Aitken’s Δ² process. This method is particularly effective for functions where the derivative is difficult to compute or is not readily available. (Wikipedia: Steffensen’s Method)

Search methods, such as BracketSearchUp and BracketSearchDown, are designed to incrementally expand or subdivide the search bounds to find a bracket where a root exists. These methods are useful when the initial guess is not near the actual root.

The following search methods are available in the library:

  • BracketSearchUp: A specialized search algorithm designed to incrementally expand the search bounds upwards (increasing values) to find a bracket where a root exists.

  • BracketSearchDown: A specialized search algorithm designed to incrementally expand the search bounds downwards (decreasing values) to find a bracket where a root exists.

  • BracketExpandUp: A specialized search algorithm designed to incrementally expand the upper bound upwards (increasing values) while keeping the lower bound fixed. This is useful for finding a bracket where a root exists when the initial guess is lower than the actual root.

  • BracketExpandDown: A specialized search algorithm designed to incrementally expand the lower bound downwards (decreasing values) while keeping the upper bound fixed. This is useful for finding a bracket where a root exists when the initial guess is higher than the actual root.

  • BracketExpandOut: A specialized search algorithm designed to incrementally expand both the lower and upper bounds symmetrically outwards. This is useful for finding a bracket where a root exists when the initial guess is not near the actual root.

  • BracketSubdivide: A specialized search algorithm designed to subdivide the current search bounds into smaller segments in an attempt to find a bracket where a root exists.

Using these methods should be done through the fsolve, fdfsolve, and search functions, which are designed to be easy to use and provide a consistent interface for all the solvers.

Important Considerations

It’s imperative to recognize that root-finding functions are designed to locate a single root at any given instance. In cases where multiple roots are present within the search range, the function will identify and return the initial root it encounters. Pinpointing which root will be found in a region with several roots is generally unpredictable. Notably, attempting to find a root in such areas usually doesn’t trigger any error messages, despite the inherent challenges.

Quick Start

The library provides the following functions for finding roots in one-dimensional functions:

  • fsolve: A function for finding a root of a one-dimensional function using a bracketing method. It takes a bracketing method object and an initial bracket around the root to be found.

  • fdfsolve: A function for finding a root of a one-dimensional function using a polishing method. It takes a polishing method object and an initial guess for the root.

  • search: A function for finding a bracket where a root exists when the initial guess is not near the actual root. It takes a search method object and an initial guess for the root.

These functions return a proxy object that contains the result of the root-finding process. The result can be accessed using the result member function, which returns the estimated root.

Note

The proxy objects returned by the root-finding functions are not intended to be stored or passed around. They are designed to be used immediately to access the result of the root-finding process. Consequently, they cannot be copied or moved, and the member functions are only available for r-value references.

The following example shows how to find the root of the function \(f(x) = x^2 - 5\) (defined as a lambda) using a bracketing method:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}).result();

Similarly, the same function can be solved using Newton’s method (with the derivative being \(f'(x) = 2x\))

using nxx::roots::fdfsolve;
using nxx::roots::Newton;
auto result = fdfsolve<Newton>([](double x){return x * x - 5;},
                               [](double x){return 2 * x;}), 1.25).result();

Note

The derivative is not required to be provided, as it will be computed numerically using a finite difference approximation. However, this can be less accurate and less efficient than providing the derivative directly.

Finally, the following example shows how to find a bracket where a root exists when the initial guess is not near the actual root:

using nxx::roots::search;
using nxx::roots::BracketSearchUp;
auto result = search<BracketSearchUp>([](double x){return x * x - 5;}, { 0.0, 1.0 }).result();

In addition to the arguments provided in the examples above, the root-finding functions also take two optional arguments for terminating the iterations, namely the relative error (epsilon) and the maximum number of iterations. Those arguments can be provided in any order; the only requirement is that epsilon must be a floating-point number and the maximum number of iterations must be an integer.

As can be seen from the previous examples, the root-finding solvers can be invoked in a single line of code, which makes them very easy to use.

Advanced Usage

Although using the root-finding functions and algorithms can be as simple as shown in the previous examples, there are a number of advanced features and considerations that can be useful in more complex scenarios.

Customizing the Solvers

In many cases, the default settings of the solvers will be sufficient. However, in some cases it may be required to provide more fine-grained control over the solvers and the criteria for terminating the iterations. This can be done by creating a “stop token” object and passing it to the solver.

A stop token is an object that provides a way to stop the iterations of the solver based on a set of criteria. It can be any callable object, such as a lambda, that takes the solver’s state as an argument and returns a boolean value indicating whether the iterations should continue or stop. The stop token can be used to check for convergence, reach a maximum number of iterations, or any other condition that should stop the iterations.

The following example shows how to create a custom stop token and use it with the fsolve function:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto stop_token = [](const auto& solver) { return solver.iter() >= 100; };
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}, stop_token).result();

In this example, the stop token is a lambda that checks whether the number of iterations has reached 100. If the condition is met, the iterations will stop and the result will be returned.

As an alternative to providing the stop token as a function argument, it can also be provided as a second template argument, after the solver algorithm. This can be useful when the stop token is a complex type or when it needs to be reused with multiple solvers:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto stop_token = [](const auto& solver) { return solver.iter() >= 100; };
using StopToken = decltype(stop_token);
auto result = fsolve<Bisection, StopToken>([](double x){return x * x - 5;}), {0.0, 2.5}).result();

Both syntaxes are equivalent, so which one to use is a matter of taste.

However, determining when to terminate the iterations is not the only use of the stop token object; it can also be used to provide additional information about the iterations, such as the number of iterations that have been performed. This can be useful for logging or debugging purposes:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto stop_token = [](const auto& solver) {
    std::cout << "Iteration: " << solver.iter() << " | ";
    std::cout << "Guess:     " << solver.guess() << std::endl;
    return solver.iter() >= 100;
    return false; };
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}, stop_token).result();

This example shows how the stop token can be used to print the iteration number and the current guess at each iteration. This can be useful for understanding how the solver is converging and for diagnosing any issues that may arise. The output can be redirected to a file or a log for further analysis.

Finally, a stop token can also be used to deal with errors. For example, if reaching the maximum number of iterations is considered an error, the stop token can be used to throw an exception when the condition is met. Or if the current guess of the root is NaN or infinite, the stop token can be used to throw an exception as well:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto stop_token = [](const auto& solver) {
    if (solver.iter() >= 100) {
        throw std::runtime_error("Maximum number of iterations reached");
    }
    if (std::isnan(solver.guess()) || std::isinf(solver.guess())) {
        throw std::runtime_error("Invalid root");
    }
    return false; };
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}, stop_token).result();

As the previous examples show, the stop token can be used to provide fine-grained control over the solvers and the criteria for terminating the iterations. This can be useful in a wide range of scenarios, from simple cases where the default settings are sufficient to complex cases where more control is needed.

Customizing the Output

In addition to customizing the solvers, it is also possible to customize the output of the root-finding functions. As mentioned earlier, the solver functions return a proxy object that contains the result of the root-finding process. This proxy object can be used to access the final result of the root-finding, using the result member function.

In some cases, however, it may be useful to access additional information about the result, or to get the result in a different format. In the simplest form, a type that can be constructed from a floating point number, can be passed as a template argument to the result function:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}).result<double>();

In more complex scenarios, a callable object, such as a lambda, can be passed as an argument. This can be usefule in situations where some error condition need to be reported, but where throwing an exception is not appropriate. This could be the C++23 std::expected type, or one of the alternatives such as tl::expected from TartanLlama:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
    auto outputter = [](const auto &data) -> tl::expected<decltype(data.guess), std::string> {
    auto [iter, guess, previous] = data;
    using expected = tl::expected<decltype(guess), std::string>;

    if (iter >= 5) return tl::make_unexpected(std::string("Too many iterations"));
    return expected(guess);
};
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}).result(outputter);

The outputter can also be provided as a template argument. The functioanlity is the same, so which one to use is a matter of taste:

using nxx::roots::fsolve;
using nxx::roots::Bisection;
using Outputter = decltype(outputter);
auto result = fsolve<Bisection>([](double x){return x * x - 5;}), {0.0, 2.5}).result<Outputter>();

The outputter can be used to provide additional information about the result, such as the number of iterations that were performed, or the reason for stopping the iterations. This can be useful for logging or debugging purposes, or for providing more context about the result of the root-finding process.

API Reference

This section provides a detailed reference for the root-finding functions. The reference includes the function signature, the template parameters, the function arguments, and the return type.

For each function, the reference also includes a list of the available algorithms, along with a brief description of each algorithm and its use cases.

Note

The algorithm classes are not intended to be used directly, but are provided for reference. When used in the root-finding functions, the algorithm template arguments will be deduced automatically, and should not be provided explicitly.

fsolve

The fsolve functions are used to find a root of a one-dimensional function using a bracketing method. It takes a bracketing method object and an initial bracket around the root to be found. The concrete algorithm to be used is provided as a template argument, and the function to be solved is provided as a function argument. The function returns a proxy object that contains the result of the root-finding process.

Four overloads of the fsolve function are provided, with different combinations of arguments, enabling the user to provide custom termination conditions and output formats.

template<template<typename, typename> class SOLVER_T, IsFloatOrComplexInvocable FN_T, IsFloatStruct STRUCT_T, typename ...ARGS>
auto nxx::roots::fsolve(FN_T func, STRUCT_T bounds, ARGS... args)

A function template that solves a root-finding problem using a specified solver and optionally the required epsilon value and/or maximum iteration count.

This function template is an overload of the fsolve function that accepts a structure for the bounds, such as a std::pair, a std::tuple, or a custom structure with two floating point elements. It uses a specified solver to find the root of a given function within the provided bounds. The required epsilon value and/or maximum iteration count can be optionally passed as additional arguments. If no additional arguments are passed, the default values for the epsilon value and maximum iteration count are used.

Template Parameters:

  • SOLVER_T – The type of the solver to use for the root-finding problem. Must be one of the pre-defined solver classes, such as Bisection, Ridder, or RegulaFalsi, or a custom solver class that meets the requirements.

  • FN_T – The type of the function for which the root is being found. Must be invocable with a floating point or complex number. This can be a function pointer, a function object, or a lambda function.

  • STRUCT_T – The type of the structure for the bounds. This can be a std::pair, std::tuple, or a custom structure. The structure must have exactly two floating point elements.

  • ARGS – The type of the additional arguments. These arguments are passed to the fsolve_common function.

Parameters:

  • func – The function for which the root is being found.

  • bounds – The bounds within which the root is being found. Must be a structure, such as a std::pair or a std::tuple, that contains exactly two floating point elements.

  • args – The additional arguments passed to the fsolve_common function.

Returns:

The result of the root-finding problem.

Pre:

The function must be continuous and have opposite signs at the bounds to ensure that a root exists within the bounds.

Post:

The result of the root-finding problem is returned as a BracketSolverResult object, which can be used to retrieve the result in a specified format.

template<template<typename, typename> class SOLVER_T, typename TOKEN_T, IsFloatOrComplexInvocable FN_T, IsFloatStruct STRUCT_T>
auto nxx::roots::fsolve(FN_T func, STRUCT_T bounds)

A function template that solves a root-finding problem using a specified solver and termination condition.

This function template is an overload of the fsolve function that accepts a structure for the bounds, such as a std::pair, a std::tuple, or a custom structure with two floating point elements. It uses a specified solver to find the root of a given function within the provided bounds. The termination condition for the solver is specified by a callable object, that is called for each iteration of the solver. The termination condition callable must accept an IterData object, which contains the current iteration count and the current bounds around the root.

Note

This function requires that the termination condition callable is invocable with an IterData object.

Template Parameters:

  • SOLVER_T – The type of the solver to use for the root-finding problem. Must be a template class that accepts two type parameters.

  • TOKEN_T – The type of the callable object that specifies the termination condition for the solver.

  • FN_T – The type of the function for which the root is being found. Must be invocable with a floating point or complex number.

  • STRUCT_T – The type of the structure for the bounds.

Parameters:

  • func – The function for which the root is being found.

  • bounds – The bounds within which the root is being found.

Returns:

A proxy object containing the result of the root-finding problem.

template<template<typename, typename> class SOLVER_T, IsFloatOrComplexInvocable FN_T, IsFloat ARG_T, size_t N, typename ...ARGS>
auto nxx::roots::fsolve(FN_T func, const ARG_T (&bounds)[N], ARGS... args)

A function template that solves a root-finding problem using a specified solver and termination condition.

This function template is an overload of the fsolve function that accepts a fixed-size array for the bounds. It uses a specified solver to find the root of a given function within the provided bounds. The termination condition for the solver is specified by a callable object.

Note

This function requires that the termination condition callable is invocable with an IterData object.

Template Parameters:

  • SOLVER_T – The type of the solver to use for the root-finding problem. Must be a template class that accepts two type parameters.

  • FN_T – The type of the function for which the root is being found. Must be invocable with a floating point or complex number.

  • ARG_T – The type of the argument to the function. Must be a floating point number.

  • N – The size of the array for the bounds. Must be 2.

  • ARGS... – The type of the additional arguments. These arguments are passed to the fsolve_common function.

Parameters:

  • func – The function for which the root is being found.

  • bounds – The bounds within which the root is being found. Must be a fixed-size array of 2 elements.

  • args – The additional arguments passed to the fsolve_common function.

Returns:

The result of the root-finding problem.

template<template<typename, typename> class SOLVER_T, typename TOKEN_T, IsFloatOrComplexInvocable FN_T, IsFloat ARG_T, size_t N>
auto nxx::roots::fsolve(FN_T func, const ARG_T (&bounds)[N])

A function template that solves a root-finding problem using a specified solver and termination condition.

This function template is an overload of the fsolve function that accepts a fixed-size array for the bounds. It uses a specified solver to find the root of a given function within the provided bounds. The termination condition for the solver is specified by a callable object.

Note

This function requires that the termination condition callable is invocable with an IterData object.

Template Parameters:

  • SOLVER_T – The type of the solver to use for the root-finding problem. Must be a template class that accepts two type parameters.

  • TOKEN_T – The type of the callable object that specifies the termination condition for the solver.

  • FN_T – The type of the function for which the root is being found. Must be invocable with a floating point or complex number.

  • ARG_T – The type of the argument to the function. Must be a floating point number.

  • N – The size of the array for the bounds. Must be 2.

Parameters:

  • func – The function for which the root is being found.

  • bounds – The bounds within which the root is being found. Must be a fixed-size array of 2 elements.

Returns:

The result of the root-finding problem.

Algorithms

The following algorithms are available for use with the fsolve functions. Note that the algorithms are not intended to be used directly, but are provided for reference. Also, when used in the fsolve function, the algorithm template arguments will be deduced automatically, and should not be provided explicitly.

template<IsFloatInvocable FN, IsFloat ARG_T = double>
class Bisection : public nxx::roots::detail::BracketingBase<Bisection<FN, double>, FN, double>

Defines the Bisection class for performing the bisection method of root bracketing.

The Bisection class template is an implementation of the classic bisection method for root finding. This method is a bracketing algorithm that repeatedly bisects an interval and then selects a subinterval in which a root must lie for further processing. It inherits from a base class that provides common functionalities for root bracketing algorithms, and adds the specific iteration logic for the bisection method. The class is templated to accept a function and an optional argument type.

Template Parameters:

  • FN – The type of the function for which the root is being bracketed.

  • ARG_T – The type of the argument to the function, defaults to double.

template<IsFloatInvocable FN, IsFloat ARG_T = double>
class Ridder : public nxx::roots::detail::BracketingBase<Ridder<FN, double>, FN, double>

Defines the Ridder class for performing Ridder’s method of root bracketing.

Ridder’s method is a root-finding algorithm that provides a more robust and often faster convergence than simple bisection. This class template, Ridder, inherits from a base class that provides common functionalities for root bracketing algorithms, and adds the specific iteration logic for Ridder’s method. It is templated to accept a function and an optional argument type.

Template Parameters:

  • FN – The type of the function for which the root is being bracketed.

  • ARG_T – The type of the argument to the function, defaults to double.

template<IsFloatInvocable FN, IsFloat ARG_T = double>
class RegulaFalsi : public nxx::roots::detail::BracketingBase<RegulaFalsi<FN, double>, FN, double>

Defines the RegulaFalsi class for performing the regula falsi (false position) method of root bracketing.

The RegulaFalsi class template implements the regula falsi method, also known as the false position method, for root finding. This method is a bracketing algorithm similar to the bisection method but, instead of bisecting the interval, it uses a linear approximation to guess the root. This class inherits from a base class that provides common functionalities for root bracketing algorithms and adds the specific iteration logic for the regula falsi method. It is templated to accept a function and an optional argument type.

Template Parameters:

  • FN – The type of the function for which the root is being bracketed.

  • ARG_T – The type of the argument to the function, defaults to double.

fdfsolve

template<template<typename, typename, typename> class SOLVER_T, IsFloatOrComplexInvocable FN_T, IsFloatOrComplexInvocable DERIV_T, IsFloatOrComplex GUESS_T, typename ...ARGS>
auto nxx::roots::fdfsolve(FN_T func, DERIV_T derivative, GUESS_T guess, ARGS... args)

A function template that solves a root finding problem using a specified solver.

This function template, fdfsolve, provides a generic implementation for root finding algorithms that utilize polishing solvers. It is designed to work with solvers that conform to the requirements of polishing solvers, such as having a defined IsPolishingSolver static member, initialization, and iteration methods. The function handles initialization, iteration, and convergence checking, returning the result along with any potential errors encountered during the solving process.

Template Parameters:

  • SOLVER_T – The template class of the solver to be used. Must be a valid polishing solver type.

  • FN_T – The type of the function for which the root is being refined.

  • DERIV_T – The type of the derivative function of FN_T.

  • GUESS_T – The type of the initial guess for the root.

  • ARGS – The type of additional arguments passed to the function.

Parameters:

  • func – The function for which the root is being refined.

  • derivative – The derivative of the function.

  • guess – The initial guess for the root.

  • args – Additional arguments passed to the function.

Returns:

The result of the root finding process.

template<template<typename, typename, typename> class SOLVER_T, IsFloatOrComplexInvocable FN_T, IsFloatOrComplex GUESS_T, typename ...ARGS>
auto nxx::roots::fdfsolve(FN_T func, GUESS_T guess, ARGS... args)
template<template<typename, typename, typename> class SOLVER_T, typename TOKEN_T, IsFloatOrComplexInvocable FN_T, IsFloatOrComplexInvocable DERIV_T, IsFloatOrComplex GUESS_T>
auto nxx::roots::fdfsolve(FN_T func, DERIV_T derivative, GUESS_T guess)

A function template that solves a root finding problem using a specified solver and termination token.

This function template, fdfsolve, provides a generic implementation for root finding algorithms that utilize polishing solvers. It is designed to work with solvers that conform to the requirements of polishing solvers, such as having a defined IsPolishingSolver static member, initialization, and iteration methods. The function handles initialization, iteration, and convergence checking, returning the result along with any potential errors encountered during the solving process.

Note

This function requires the termination token to be a callable object that accepts an IterData object.

Template Parameters:

  • SOLVER_T – The template class of the solver to be used. Must be a valid polishing solver type.

  • TOKEN_T – The type of the termination token. Must be a callable object that accepts an IterData object.

  • FN_T – The type of the function for which the root is being refined.

  • DERIV_T – The type of the derivative function of FN_T.

  • GUESS_T – The type of the initial guess for the root.

Parameters:

  • func – The function for which the root is being refined.

  • derivative – The derivative of the function.

  • guess – The initial guess for the root.

Returns:

The result of the root finding process.

template<template<typename, typename, typename> class SOLVER_T, typename TOKEN_T, IsFloatOrComplexInvocable FN_T, IsFloatOrComplex GUESS_T>
auto nxx::roots::fdfsolve(FN_T func, GUESS_T guess)

Algorithms

The following algorithms are available for use with the fdfsolve functions:

template<IsFloatOrComplexInvocable FN, IsFloatOrComplexInvocable DFN, IsFloatOrComplex ARG_T = double>
class Newton : public nxx::roots::detail::PolishingBase<Newton<FN, DFN, double>, FN, DFN, double>

Defines the Newton class for performing Newton’s method root polishing.

The Newton class template is a specialized implementation of the root polishing algorithm known as Newton’s method (or the Newton-Raphson method). It inherits from a base class that provides common functionalities for root polishing algorithms, and adds the specific iteration logic for Newton’s method. This class is templated to accept a function, its derivative, and an optional argument type.

Template Parameters:

  • FN – Type of the function for which the root is being refined.

  • DFN – Type of the derivative function of FN.

  • ARG_T – Type of the argument to the function and its derivative, defaults to double.

template<IsFloatOrComplexInvocable FN, IsFloatOrComplexInvocable DFN, IsFloatOrComplex ARG_T = double>
class Secant : public nxx::roots::detail::PolishingBase<Secant<FN, DFN, double>, FN, DFN, double>

Defines the Secant class for performing the secant method of root polishing.

The Secant class template implements the secant method, an iterative root finding algorithm that uses a succession of roots of secant lines to approximate a root of a function. This class template inherits from a base class that provides common functionalities for root polishing algorithms, and adds specific logic for the secant method iterations. It is designed to handle an initial guess and a subsequent series of approximations to converge on the root. The secant method is a derivative-free alternative to Newton’s method and is particularly useful when the derivative of the function is not readily available or is expensive to compute.

Template Parameters:

  • FN – The type of the function for which the root is being polished.

  • DFN – The type of the derivative function of FN, used in the initial step.

  • ARG_T – The type of the argument to the function, defaults to double.

template<IsFloatOrComplexInvocable FN, IsFloatOrComplexInvocable DFN, IsFloatOrComplex ARG_T = double>
class Steffensen : public nxx::roots::detail::PolishingBase<Steffensen<FN, DFN, double>, FN, DFN, double>

Defines the Steffensen class for performing Steffensen’s method of root polishing.

Steffensen’s method is an iterative root finding algorithm that improves upon the simple fixed-point iteration by incorporating a form of Aitken’s Δ² process. This class template, Steffensen, inherits from a base class that provides common functionalities for root polishing algorithms, and adds the specific logic for Steffensen’s method iterations. It starts with a Newton-Raphson step and then switches to Steffensen’s method for subsequent iterations. This method is particularly effective for functions where the derivative is difficult to compute or is not readily available.

Template Parameters:

  • FN – The type of the function for which the root is being polished.

  • DFN – The type of the derivative function of FN, used in the initial step.

  • ARG_T – The type of the argument to the function, defaults to double.

Design and Implementation Details

The solver functions documented above should be sufficient for the vast majority of use cases. However, in some cases it may be necessary to use the solvers directly, or to create custom solvers. This section provides an overview of the design and implementation details of the solvers, and how they can be used and customized.

Interface

Both the bracketing and polishing algorithms are implemented using the overall architecture: a base class is defined, for keeping track of internal state between iterations, and the individual algorithms inherits from the base class. However, in order to avoid virtual functions, the architecture is implemented using static polymorphism through the Curiously Recurring Template Pattern (CRTP) [1].

The purpose of this design is that it avoids dynamic polymorphim and virtual functions, while ensuring that the solvers share a common interface. The downside of this approach is that all dependencies has to be resolved during compile-time, and it is not possible to dynamically plug in a different solver.

While the base classes (BracketingBase and PolishingBase, respectively) will not be called directly in client code, it is useful to know what the classes look like, as the individual solvers will inherit the interface of the base classes.

BracketingBase

The BracketingBase class (located in the root::detail namespace) look as follows:

template<typename DERIVED, typename FUNCTION_T, typename ARG_T>
class BracketingBase

Provides a base class template for root bracketing algorithms.

The BracketingBase class template serves as a foundational component for algorithms that bracket roots of a given function. It encapsulates common functionalities such as storing the objective function and maintaining the current bounds around the root. This class enforces certain type constraints on the template parameters to ensure compatibility with root bracketing algorithms.

Template Parameters:

  • DERIVED – The subclass inheriting from BracketingBase.

  • FUNCTION_T – The type of the function for which the root is being bracketed.

  • ARG_T – The type of the argument to the function.

Subclassed by nxx::roots::Bisection< FN, ARG_T >, nxx::roots::RegulaFalsi< FN, ARG_T >, nxx::roots::Ridder< FN, ARG_T >

Public Functions

inline RESULT_T evaluate(IsFloat auto value)

Evaluates the function at a given value.

Parameters:

value – The value at which the function is to be evaluated.

Returns:

The result of evaluating the function at the specified value.

Private Members

FUNCTION_T m_func = {}

The function object to find the root for.

BOUNDS_T m_bounds = {}

Holds the current bounds around the root.

PolishingBase

Similarly, the PolishingBase class (located in the root::detail namespace) look as follows:

template<typename SUBCLASS, typename FUNCTION_T, typename DERIV_T, typename ARG_T>
class PolishingBase

Provides a base class template for root polishing algorithms.

The PolishingBase class template serves as a foundational component for algorithms that refine or ‘polish’ roots of a given function. It encapsulates common functionalities such as storing the objective function, its derivative, and the current guess of the root. This class enforces certain type constraints on the template parameters to ensure compatibility with root polishing algorithms.

Template Parameters:

  • SUBCLASS – The subclass inheriting from PolishingBase.

  • FUNCTION_T – The type of the function for which the root is being polished.

  • DERIV_T – The type of the derivative function of FUNCTION_T.

  • ARG_T – The type of the argument to the function and its derivative.

Subclassed by nxx::roots::Newton< FN, DFN, ARG_T >, nxx::roots::Secant< FN, DFN, ARG_T >, nxx::roots::Steffensen< FN, DFN, ARG_T >

Public Functions

template<typename T>
inline auto evaluate(T value)

Evaluates the function with the given value.

Template Parameters:

T – Type of the value, must be a float or complex type.

Parameters:

value – The value to evaluate the function at.

Returns:

The result of the function evaluation.

template<typename T>
inline auto derivative(T value)

Evaluates the derivative function with the given value.

Template Parameters:

T – Type of the value, must be a float or complex type.

Parameters:

value – The value to evaluate the derivative at.

Returns:

The result of the derivative function evaluation.

Private Members

FUNCTION_T m_func = {}

The function object to find the root for.

DERIV_T m_deriv = {}

The function object for the derivative.

RESULT_T m_guess

The current root estimate.

Implementation

The implementation of the solvers is based on the algorithms described in the previous sections. The solvers are implemented as classes, with the individual algorithms being implemented through the function call operator (operator()). This makes the implementation reasonably straightforward and easy to understand.