Roots.jl

Root finding methods are characterized by their order of convergence $q$ (where $e_{i+1} \approx e_i^q$) and their asymptotic efficiency $q^{1/p}$, where $p$ is the number of function evaluations per step.

When choosing a method, consider the trade-off between the convergence rate and the number of function evaluations (F evals) required per iteration. For example, Newton's method has quadratic convergence ($q=2$) but requires 2 function evaluations, whereas the Secant method has an order of $\varphi \approx 1.618$ but only requires 1 evaluation per step.

Common method types include:

• Classical: Methods like Newton, Halley, and SuperHalley that typically use derivatives.
• Derivative Free: Methods like Secant, Steffensen, and Order16 that do not require derivatives.
• Bracketing: Methods like Brent, ITP, and Ridders that work within a known interval $[a, b]$.
• Robust: Methods like Schroder or Thukral designed for stability.

| Type            | Method                       | Order                  | F evals | Asymptotic efficiency                 |
|:--------------- | :--------------------------- | :--------------------- | :------ | :------------------------------------ |
| Hybrid          | Order0                       |                        |         | ``\approx 1.618\dots``               |
| Derivative Free | Secant                       | ``\varphi=1.618\dots`` | ``1``   | ``1.618\dots``                       |
| Derivative Free | Steffensen                   | ``2``                  | ``2``   | ``1.414\dots``                       |
| Derivative Free | Order5                       | ``5``                  | ``4``   | ``1.495\dots``                       |
| Derivative Free | Order8                       | ``8``                  | ``4``   | ``1.681\dots``                       |
| Derivative Free | Order16                      | ``16``                 | ``5``   | ``1.718\dots``                       |
| Classical       | Newton                       | ``2``                  | ``2``   | ``1.414\dots``                       |
| Classical       | Halley                       | ``3``                  | ``3``   | ``1.442\dots``                       |
| Classical       | QuadraticInverse             | ``3``                  | ``3``   | ``1.442\dots``                       |
| Classical       | ChebyshevLike                 | ``3``                  | ``3``   | ``1.442\dots``                       |
| Classical       | SuperHalley                  | ``3``                  | ``3``   | ``1.442\dots``                       |
| MultiStep       | LithBoonkkampIJzerman{S,D}    | ``p^s=\sum p^k(d+\sigma_k)`` | ``D+1`` | varies, ``1.92\dots`` max            |
| Bracketing      | BisectionExact               | ``1``                  | ``1``   | ``1``                                |
| Bracketing      | A42                          | ``(2 + 7^{1/2})``      | ``3,4`` |``(2 + 7^{1/2})^{1/3} = 1.6686\dots`` |
| Bracketing      | AlefeldPotraShi              |                        | ``3,4`` | ``1.618\dots``                       |
| Bracketing      | Brent                        | ``\leq 1.89\dots``     | ``1``   | ``\leq 1.89\dots``                  |
| Bracketing      | ITP                          | ``\leq \varphi``       | ``1``   | ``\leq \varphi``                    |
| Bracketing      | Ridders                      | ``1.83\dots``           | ``2``   | ``1.225\dots``                      |
| Bracketing      | RegularFalsi{:classic}        | ``1``                  | ``1``   | ``1``                                |
| Bracketing      | RegularFalsi{:Illinois}       | ``1.442\dots``         | ``1``   | ``1.442\dots``                       |
| Bracketing      | RegulaFalsi{:AndersonBjork}  | ``1.681\dots``         | ``1``   | ``1.681\dots``                       |
| Bracketing      | RegulaFalsi{:Ford4}          | ``1.681\dots``         | ``1``   | ``1.681\dots``                       |
| Bracketing      | ModAB                        | ``≈1.7\dots``          | ``1``   | ``1.7\dots``                         |
| Bracketing      | LithBoonkkampIJzermanBracket | ``2.91``               | ``3``   | ``1.427\dots``                       |
| Robust          | King                         | ``\varphi=1.618\dots`` | ``2``   | ``1.272\dots``                       |
| Robust          | Esser                        | ``2``                  | ``3``   | ``1.259\dots``                       |
| Robust          | Schroder                     | ``2``                  | ``3``   | ``1.259\dots``                       |
| Robust          | Thukral3                     | ``3``                  | ``4``   | ``1.316\dots``                       |
| Robust          | Thukral4                     | ``4``                  | ``5``   | ``1.319\dots``                       |
| Robust          | Thukral5                     | ``5``                  | ``6``   | ``1.307\dots``                       |

When using Float64, mathematical zeros may not behave as expected due to precision limits:

• Multiple Zeros: A single mathematical zero might appear as multiple floating-point zeros due to evaluation differences (e.g., iszero(f(x)) being true for multiple adjacent values).
• Sign Changes: A mathematical function that changes sign might not satisfy f(x) * f(nextfloat(x)) < 0 exactly due to rounding. An exact zero is identified if iszero(f(x)) is true OR if there is a sign change between adjacent floating-point numbers.
• Residual Magnitude: At a floating-point zero $x$, the function value $f(x)$ is not necessarily $0$. It is typically on the scale of $f'(x) \cdot |x| \cdot \epsilon$, where $\epsilon$ is the machine epsilon.

Example: Verifying sign changes near a zero

f(x) = exp(x) - x^4;
F(x) = sign(f(x));
x = 8.613169456441398
F(prevfloat(x)), F(x), F(nextfloat(x))
# (-1.0, -1.0, 1.0)

Root finding involves balancing absolute and relative tolerances. Because errors in function evaluation can be proportional to the size of $x$, a robust check often uses both:

abs(f(x)) < max(atol, abs(x) * rtol)

Relative Tolerance Pitfalls

For functions with sublinear growth, extremely large values of $x$ might be incorrectly identified as zeros because the relative tolerance $|x| \cdot \epsilon$ becomes larger than the function value $f(x)$.

Example of misidentified zero with Thukral8:

find_zero(cbrt, 1, Roots.Thukral8())
# 1.725042287244107e23

To avoid this, you can set the relative tolerance to $0$ and use a more generous absolute tolerance (e.g., sqrt(eps()) or 1e-8).

Roots.jl Tolerance Strategy

Roots.jl uses a dual approach for faster algorithms:

1. Tight Check: Checks if the difference between the last two $x_n$ values is small AND the residual $f(x_n)$ is small with a tight tolerance.
2. Relaxed Check: If the $x$ values are close AND the function value is close to zero with a relaxed tolerance, an approximate zero is declared.

Following the pattern used in DifferentialEquations.jl, Roots provides a structured way to solve problems using ZeroProblem and the CommonSolve.jl interface. This is useful for separating the problem definition from the solver choice.

1. Define the problem: Create a ZeroProblem(f, x0) where f is the function and x0 is the initial guess or bracket.
2. Solve the problem: Use solve(Z, M) where Z is the problem and M is the algorithm (e.g., Secant(), Bisection(), Order2()).

This interface supports parameterized functions by passing parameters as positional arguments or keyword arguments to solve.

using Roots

f(x) = sin(x)
x0 = (3, 4)
M = Secant()

# 1. Setup the problem
Z = ZeroProblem(f, x0)

# 2. Solve with a specific algorithm
result = solve(Z, M)

# Solving parameterized functions
g(x, p=1) = cos(x) - x/p
Z_param = ZeroProblem(g, (0.0, pi/2))
solve(Z_param, Secant(), 2) # p=2 as positional argument
solve(Z_param, Bisection(); p=3, xatol=1/16) # p=3 via keyword