# Other zero-finding algorithms

{{< include ../_common_code.qmd >}}

This section uses these add-on packages:


```{julia}
using CalculusWithJulia
using Plots
plotly()
using Roots
using SymPy
```

---


There are *numerous* zero-finding methods in addition to the secant method and Newton's method. This section shows a few different directions. It then discusses the topic of when to terminate an algorithm.

This section is entirely optional, none of the algorithms discussed below are utilized in the sequel.

## Other methods

We discuss variations of both Newton's method and the bisection method.

### Estimating the derivative

[Sidi](https://www.math.nthu.edu.tw/~amen/2008/070227-1.pdf) starts with Newton's method with its update step

$$
x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}
$$

and notes that the secant method just uses the slope of the secant line between $x_i$ and $x_{i-1}$ to estimate $f'(x_i)$. The secant line is the *linear* polynomial interpolating the two points $(x_i, f(x_i))$ and $(x_{i-1}, f(x_{i-1}))$ and the slope the derivative of this polynomial. Sidi generalizes this to approximate the function $f'(x_i)$ by using more points from the algorithm to interpolate a polynomial at $x_i, x_{i-1}, \dots, x_{i-k}$ and then using the derivative of this polynomial to estimate the derivative of $f$ at $x_i$.


Let's consider a case with three points or $k=2$. Here we use  formula (10) of Sidi computing the derivative of the interpolating polynomial at $x_n$ specialized for $k=2$.

```{julia}
function pprime_n2(xn_2, xn_1, xn, yn_2, yn_1, yn)
    m = (yn - yn_1) / (xn - xn_1)
    m_1 = (yn_1 - yn_2) / (xn_1 - xn_2)
    m + (m - m_1)/(xn - xn_2) * (xn - xn_1)
end
```

Typically, we would start with two initial points and apply the secant method to produce a third, but for this example we will use an initial three points.

```{julia}
f(x) = cos(x) - x/2
xs = [0, pi/2, pi/4]
ys = f.(xs)                   # y₀, y₁, y₂

xx = xs[end] - ys[end] / pprime_n2(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx); push!(ys, yy)  # y₃
```

This method generalizes the secant method with a convergence rate of $1.83928\cdots$. For this problem we see it takes  5 iterations to converge to machine tolerance:

```{julia}
xx = xs[end] - ys[end] / pprime_n2(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx); push!(ys, yy) # y₄ = 0.0004182051168989398

xx = xs[end] - ys[end] / pprime_n2(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx); push!(ys, yy) # y₅ = -3.6161489780361933e-7

xx = xs[end] - ys[end] / pprime_n2(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx); push!(ys, yy) # y₆ = 3.609335053056384e-13

xx = xs[end] - ys[end] / pprime_n2(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx); push!(ys, yy) # y₇ = 0.0

xx
```

### Estimating the inverse function

Suppose $f^{-1}$ exists in a neighborhood of $\alpha$ and we have generated steps in our algorithm $x_0, x_1, \dots, x_n$. We can find $\alpha$ from $f^{-1}(0)$. Typically though, we wouldn't have the inverse function, but we can use facts about functions such as linearization, which *near* $0$ has:

$$
f^{-1}(y) \approx f^{-1}(0) + (f^{-1})'(0)\cdot(y) = \alpha + (f^{-1})'(0)\cdot y
$$

Solving for $\alpha$ gives for $(x_i, f(x_i))$

$$
\alpha \approx f^{-1}(f(x_i)) - f^{-1}(0) f(x_i).
$$

Replacing $f^{-1}(0) \approx (f^{-1})'(y_i) = (f^{-1})'(f(x_i)) = 1/f'(x_i)$ we get:

$$
\alpha \approx x_i - \frac{f(x_i)}{f'(x_i)}.
$$

Which is basically Newton's method.

The above uses one point, $(x_i, f(x_i))$ and the fact that $f$ is differentiable. What if two (or more) points were used, would that give some insight?

Here is some code that interpolates a *polynomial* for $f^{-1}(y)$ through $k$ points and then solves for it's value at $0$, which is $a_0$.


```{julia}
@syms x[1:5] y[1:5] a[0:5]
a₀ = first(a)
function I(k)
    eqs = Tuple(sum(a[j] * y[i]^(j-1) for j in 1:k) ~ x[i] for i in 1:k)
    sols = solve(eqs, Tuple(a[1:k]))
    sols[a₀]           # intercept
end
```


We use $k=2$ and see what comes:

```{julia}
I(2)
```

We can see this is a rewriting of the secant method through:

```{julia}
x1, x2 = x[1:2]; y1, y2 = y[1:2]
m = (y2 - y1) / (x2 - x1)
secant_method = x2 - (1/m) * y2
simplify(I(2) - secant_method)
```
An inverse quadratic step ($k=2$) is utilized by Brent's method, as possible, to yield a rapidly convergent bracketing algorithm implemented as a default zero finder in many software languages.  `Julia`'s `Roots` package implements the method in `Roots.Brent()`.

To see an example of inverse quadratic, we first make a function to compute the next $x$ value, given three previous $x$ and $f(x)$ values.

```{julia}
u = lambdify(I(3), (x[1:3]..., y[1:3]...))
```

Let's try initial values $(x_0, x_1, x_2) = (0, \pi/2, \pi/4)$:

```{julia}
f(x) = cos(x) - x/2
xs = [0, pi/2, pi/4]
ys = f.(xs)
```

Now we do a step. The new values is "pushed" to the vector of values.

```{julia}
xx = u(xs[end-2:end]..., ys[end-2:end]...)
yy = f(xx)
push!(xs, xx)
push!(ys, yy)
xx, yy
```

We know do a few more steps, the value of `yy` is shown as a comment.

```{julia}
xx = u(xs[end-2:end]..., ys[end-2:end]...); yy = f(xx)
push!(xs, xx); push!(ys, yy) # yy = 0.0011053827937966831

xx = u(xs[end-2:end]..., ys[end-2:end]...); yy = f(xx)
push!(xs, xx); push!(ys, yy) # yy = -2.113895167021873e-6

xx = u(xs[end-2:end]..., ys[end-2:end]...); yy = f(xx)
push!(xs, xx); push!(ys, yy) # yy = 1.1904810470753091e-11

xx = u(xs[end-2:end]..., ys[end-2:end]...); yy = f(xx)
push!(xs, xx); push!(ys, yy) # yy = 0.0

xx
```

The convergence happens quite quickly with this well-behaved problem.

An inverse cubic interpolation is utilized by [Alefeld, Potra, and Shi](https://dl.acm.org/doi/10.1145/210089.210111) which gives an asymptotically even more rapidly convergent algorithm than Brent's (implemented in `Roots.AlefeldPotraShi()` and also `Roots.A42()`). This is used as a finishing step in many cases by the default hybrid `Order0()` method of `find_zero`.




### Steffensen's method

Another alternative to the secant method is Steffensen's method.

The secant method has super-linear convergence, but not quadratic convergence. It uses these points to evaluate the values $(x_i, f(x_i)$ and $(x_{i-1}, f(x_{i-1}))$. When $x_i$ converges to $\alpha$, $x_i - x_{i-1}$ will converge to $0$. The secant lines used are eventually "converging" to tangent lines.

Steffensen's method takes a different pair of points to use for a secant line, these being $(x_n,f(x_n))$ and  $(x_n + f(x_n), f(x_n + f(x_n)))$. When $x_i \rightarrow \alpha$ it follows for a continuous $f(x)$ that $f(x_i) \rightarrow 0$, so the secant lines used by Steffensen's method will also be close to the tangent line.

[This note](https://fractal.math.unr.edu/~ejolson/701-12/code/hw2sol/hw2sol.pdf) shows that with $\eta_i$ and $\xi_i$ being values that *converge* to $\alpha$, that

$$
e_{i+1} = -e_{i}^2 \cdot
\left(\frac{f''(\xi_i)\left(f'(x_i) - \frac{1}{2} f''(\eta_i) e_i\right) +f''(\eta_i)}{
            2f'(x_i) + f''(\xi_i)f(x_i)}\right)
$$

As the following ratio converges to something non zero under assumptions, the Steffensen method has quadratic convergence.

$$
e_{i+1}/e_i^2 \rightarrow \frac{f''(\alpha)(1 + f'(\alpha))}{2f'(\alpha)}
$$


Like Newton's method this method requires $2$ function evaluations per step, but unlike Newton's method is derivative free. Steffensen's is implemented in the `Roots` package through `Roots.Steffensen()`. Steffensen's method is more sensitive to the initial guess than other methods, so in practice must be used with care, though it is a starting point for many higher-order derivative-free methods.



### Alternative bracketing methods

The bisection method has several advantages, primarily it is guaranteed to converge regardless of any assumptions on the shape of the function. It's implementation in `Roots` can handle any $x$ values as long as the function value has a sign (not `NaN` and not an error). However, it is slow---linearly convergent. There can be improvements.

#### Regula falsi

One alternative is the (modified) *regula falsi* method which replaces the midpoint ($x_i/2 + x_{i-1}/2$) with the intersection point of the line between two bracketing points $(x_i, f(x_i))$ and $(x_{i-1}, f(x_{i-1}))$ given by solving the following, which comes from similar triangles:

```{julia}
@syms xᵢ xᵢ₋₁ yᵢ yᵢ₋₁ x
only(solve(yᵢ / (x - xᵢ) ~ -yᵢ₋₁ / (xᵢ₋₁ - x), x))
```


As seen earlier, this formula is a single step of the secant method, but unlike the secant method, for this method the two points chosen to continue are picked to ensure $x_i, x_{i-1}$ form a bracketing interval.

Despite being related to the secant method, the convergence rate of *regula falsi* is only linear. Some function shapes preference a certain endpoint, whereas the secant method chooses the last two values.

@fig-regula-false-convex show that some function shapes result in  one end point being fixed as the algorithm progresses which can lead to linear convergence.

::: {#fig-regula-false-convex}
```{julia}
#| echo: false
let
    gr()
    Δ = 0.2
    function add_line!(xs, ys, f, a)
        x = (xs[1]*ys[2] - xs[2]*ys[1])/(ys[2] - ys[1])
        scatter!(collect(zip(xs, ys)); marker=(5, :orange))
        plot!(collect(zip(xs, ys)); line=(1, :gray50))
        scatter!([(x, 0)]; marker=(5, :blue))
        annotate!([(x, 0, text(a, :top))])
        plot!(plt, [(x, 0), (x, Δ)]; line=(1, :gray50))
        xs[2] = x
        ys[2] = f(x)
    end

    f(x) = 10*log(x)/x^3
    plt = plot(; empty_style..., xlims=(0.7, 3.2))
    plt = plot!(plt, f, 0.8, 3)
    plot!([(0.8, 0), (3.2, 0)]; line=(1, :black), arrow=true, side=:right)

    x₀ = 2.75
    x₁ = 0.85
    α = 1.0
    xs = [x₁, x₀]
    ys = f.(xs)

    add_line!(xs, ys, f, L"x_2")
    add_line!(xs, ys, f, L"x_3")
    add_line!(xs, ys, f, L"x_4")
    add_line!(xs, ys, f, L"x_5")
    add_line!(xs, ys, f, L"x_6")
    add_line!(xs, ys, f, L"x_7")

    plot!(plt, [(x₀, 0), (x₀, Δ)]; line=(1, :gray50))
    plot!(plt, [(x₁, 0), (x₁, Δ)]; line=(1, :gray50))
    annotate!(plt, [(x₀, 0, text(L"x_0", :top)),
                    (x₁, 0, text(L"x_1", :top)),
                    (α, 0, text(L"\alpha", :left, :top))])
    scatter!(plt, [(α, 0)]; marker=(5, :green))

    plotly()
    plt
end
```

Plot illustrating that the *regula falsi* method may have a fixed endpoint for some convex functions
:::


#### Modified regula falsi

@fig-modified-regula-falsi shows a scenario where the secant line between $(x_{i-1}, f(x_{i-1}))$ and $(x_i, f(x_i))$ crosses the $x$ axis at $x_{i+1}$ which is to the *right* of the zero $\alpha$, as it always will be for this curve and these points. A modified *regula falsi* method modifies the fixed end by using $\tilde{f}(x_i)$ and not $f(x_i)$ to compute the secant line, where $\tilde{f}$ is some multiple, $\gamma$, of $f$. In the figure, $\gamma$ is shown so that the *next* choice ($x_{i+2}$ would be its label) is exactly $\alpha$. And value for the multiplier less than this $\gamma$ will shift the intersection point to the other side of $\alpha$. The value of $\gamma$ is the ratio of the secant line slopes between $x_{i+1}$ and $\alpha$ and between $x_{i-1}$ and $\alpha$. Some choices for $\gamma$ lead to super-linear convergence.

::: {#fig-modified-regula-falsi}
```{julia}
#| echo: false
let
    gr()
    dd(f, a, b) = (f(a) - f(b)) / (a - b)
    f(x) = (x-3)^2 - 1
    plt = plot(; xlims=(0.5, 3.25), empty_style...)
    plot!(plt, [(1,0), (3.25,0)]; arrow=true, side=:right, line=(1, :gray))
    plot!(f, 1, 3.25; line=(1, :black))
    xᵢ₋₁ , xᵢ = 1.25, 2.75
    xᵢ₊₁ = (xᵢ₋₁ * f(xᵢ) -  xᵢ * f( xᵢ₋₁)) / (f(xᵢ)  - f(xᵢ₋₁))

    α = 2
    γ = dd(f, xᵢ₊₁, α) / dd(f, xᵢ₋₁, α)

    plot!(plt, [(xᵢ₋₁, 0), (xᵢ₋₁, f(xᵢ₋₁))]; line=(1, :dot))
    plot!(plt, [(xᵢ₋₁, f(xᵢ₋₁)) , (xᵢ, f(xᵢ))]; line=(1, :dash, :blue))
    plot!(plt, [(xᵢ₋₁, γ*f(xᵢ₋₁)) , (xᵢ₊₁, f(xᵢ₊₁))]; line=(1, :dash, :blue))

    scatter!(plt, [(xᵢ₋₁, γ*f(xᵢ₋₁)), (xᵢ₋₁, f(xᵢ₋₁)),
                   (α, 0), (xᵢ₊₁, 0),
                   (xᵢ₊₁, f(xᵢ₊₁)), (xᵢ, f(xᵢ))
                   ]; marker=(3, :orange))

    annotate!(plt, [
        (xᵢ₋₁, 0, text(L"x_{i-1}", :top)),
        (xᵢ₋₁, γ*f(xᵢ₋₁), text(L"\gamma \cdot f(x_{i-1})", :right)),
        (xᵢ₋₁, f(xᵢ₋₁), text(L"f(x_{i-1})", :right)),
        (α, 0, text(L"\alpha", :bottom)),
        (xᵢ₊₁, 0, text(L"x_{i+1}", :top)),
        (xᵢ, 0, text(L"x_{i}", :top))
        ])

    plotly()
    plt
end
```

Modified *regula falsi* method illustration. When midpoint $x_{i+1}$ is on same side of zero $\alpha$ as $x_i$ the *next* step will be between $x_{i-1}$ and $x_{i+1}$. *Were* $x_{i-1}$ modified by $\gamma$ the next midpoint would be an exact zero. If multiplied by a value less, then the midpoint moves to other side of $\alpha$ and would break the repeated choice of a fixed side when keeping a bracketing interval.
:::


#### Anderson Bjork

There are numerous modifications of the *regula falsi* algorithm that
employ different scaling values, we discuss one now. The
[Anderson-Bjork](https://iopscience.iop.org/article/10.1088/1757-899X/1276/1/012010/pdf)
method is a modification of the *regula falsi* method that avoids the
linear convergence when one endpoint is always fixed.

The modification works as follows, suppose the bracketing interval is $[a,b]$ and $c$ is the point found by the secant line. Then if $f(a)$ and $f(c)$ have the same sign **and** the previous step kept the right side point ($b$) fixed, then instead of using $(c, f(c))$ and $(b, f(c))$ as the new points (as $[c,b]$ is a bracket) use $(b, \gamma \cdot f(b))$ where $\gamma = 1 - f(c)/f(a)$ if $\gamma$ is positive and $\gamma=1/2$ if not. This will modify the next step in the algorithm. The $\gamma$ factors are multiplied each time, so that eventually the $y$ value used at the fixed side should lead to a midpoint on the other side of the zero, as happens when the modified value of $f(x_1)$ and $f(x_5)$ are used to find the midpoint in @fig-anderson-bjork-trajectory.

::: {#fig-anderson-bjork-trajectory}
```{julia}
#| echo: false
let
    gr()

    midpt(a,b,fa,fb) = (a*fb - b*fa)/(fb-fa)
    side = nothing
    function ABstep!(xs, ys, side, label)
        a, b = xs; fa, fb = ys
        c = midpt(a,b,fa,fb)
        fc = f(c)

        if sign(fa) == sign(fc)
            xs[1] = c
            ys[1] = fc
            if side == :right
                m = 1 - fc/fa
                m = m < 0 ? 1/2 : m
                ys[2] *= m
            else
                side = :right
            end
        else
            xs[2] = c
            ys[2] = fc
            if side == :left
                m = 1 - fc/fb
                m = m < 0 ? 1/2 : m
                ys[1] *= m
            else
                side = :left
            end
        end

        plot!([(c,0), (c,fc)]; line=(1, :dash, :gray50))
        scatter!([(c,0)]; marker=(5, :blue))
        scatter!(collect(zip(xs, ys)); marker=(5, :orange))
        plot!(collect(zip(xs, ys)); line=(1, :gray50))

        annotate!([(c, 0, text(label, :top, :left))])
        side
    end

    f(x) = 10*log(x)/x^3
    plt = plot(f, 0.8, 3; xlims=(0.7, 3.2), empty_style...)
    plot!([(0.8, 0), (3.2, 0)]; line=(1, :black), arrow=true, side=:right)

    α = 1
    x₀, x₁ = 2.75, 0.85
    xs = [x₀, x₁]
    ys = f.(xs)
    scatter!(plt, collect(zip(xs, ys)); marker=(5, :orange))
    plot!(plt, collect(zip(xs, ys)); line=(1, :gray50))
    plot!(plt, [(x₀,0), (x₀, f(x₀))]; line=(1, :dash, :gray50))
    plot!(plt, [(x₁,0), (x₁, f(x₁))]; line=(1, :dash, :gray50))


    side = ABstep!(xs, ys, side, L"x_2")
    side = ABstep!(xs, ys, side, L"x_3")
    side = ABstep!(xs, ys, side, L"x_4")
    side = ABstep!(xs, ys, side, L"x_5")
    side = ABstep!(xs, ys, side, L"x_6")

    c = midpt(xs..., ys...)
    scatter!([(c,0)];  marker=(5, :blue))
    scatter!([(α,0)];  marker=(5, :green))
    annotate!([(c,0, text(L"x_7", :top, :left)),
               (x₀, 0, text(L"x_0", :top, )),
               (x₁, 0, text(L"x_1", :top, :right)),
               (α, 0, text(L"\alpha", :top, :left))
               ])

    plotly()
    current()
end
```

Illustration of Anderson-Bjork algorithm. The point $x_1$ stays as the left-hand endpoint up until $x_6$, but by modifying $f(x_0)$ the algorithm converges super-linearly towards $\alpha$, as compared to @fig-regula-false-convex.

:::


::: {.callout-note}
## Hybrid algorithms

There are a few, newer, *hybrid algorithms* where some dynamic choice is made as to what update step should be chosen. One due to Chandrapatla (implemented in `Roots.Chandrapatla`) is a bracketing algorithm which chooses between an inverse quadratic step or a bisection step using a certain inequality. We note another, a bracketing algorithm due to Ganchovski and Traykov (with improvements by some `Julia` programmers since its inclusion in the `NonLinearSolve.jl` package) that chooses between bisection or the Anderson-Bjork update based on an estimate of how "straight" the curve is. This is implemented in `Roots.ModAB`. The latter is quite efficient over a wide range of problems.
:::




##### Examples

The function $f(x) = (x^2 + 1) \cdot \sin(x) - e^{\sqrt{\lvert x\rvert}} \cdot (x - 1) \cdot (x^2 - 5)$ has a zero *between* $0$ and $1$, and *near* $0.8$. We see how to find it using various algorithms implemented in `Roots`.

```{julia}
f(x) = (x^2 + 1) * sin(x) - exp(sqrt(abs(x))) * (x - 1) * (x^2 - 5)
x0 = 0.8
xs = (0, 1)
```

The Steffensen method needs a *nearby* estimate:

```{julia}
find_zero(f, x0, Roots.Steffensen())   # 5 iterations, 12 function evaluations
```

The `Sidi(2)` method needs a *nearby* estimate or an initial two points for a secant line which it bootstraps to get a third point. We use the bracketing interval below:

```{julia}
find_zero(f, xs, Roots.Sidi(2))        # 3 iterations, 6 function evaluations
```

For some steps, `Brent` and `Chandrapatla` use a quadratic inverse calculation, whereas `A42` uses a cubic inverse calculation:

```{julia}
find_zero(f, xs, Roots.Brent())        # 16 iterations, 18 function evaluations
find_zero(f, xs, Roots.Chandrapatla()) # 20 iterations, 22 function evaluations
find_zero(f, xs, Roots.A42())          # 4 iterations, 10 function evaluations
```

Finally, we compare *regula falsi* variants:

```{julia}
find_zero(f, xs, Roots.RegulaFalsi(:classic))       # 10 iterations, 13 function evaluations
find_zero(f, xs, Roots.RegulaFalsi(:AndersonBjork)) # 6 iterations, 9 function evaluations
find_zero(f, xs, Roots.ModAB())                     # 6 iterations, 8 function evaluations
```

For this problem, all methods converge to the same zero, but from the counts of iterations and function evaluations they differ in how the get there.


## Tolerances

Iterative zero-finding algorithms may mathematically converge, but when implemented on the computer a stopping rule must be articulated. Typically these involve the following:

* stop (and fail) if too many steps are taken
* stop when $\lvert x_i - x_{i-1} \rvert$ is quite small (as then the algorithm stops improving)
* stop when $f(x_i)$ is quite small, as it is close to being zero.

Small on the computer is a *relative* term and requires a bit of discussion.

When $\lvert x_i - x_{i-1} \rvert$ is small, we have to recall that the gap between floating point numbers depends on the size of the number, and doubles in going from $[2^{i-1}, 2^i)$ to $[2^i, 2^{i+1})$. As such, a relative tolerance is often chosen so that *small* really means that for some $\epsilon$

$$
\lvert x_i - x_{i-1} \rvert \leq \max(\lvert x_i \rvert, \lvert x_{i-1} \rvert) \cdot \epsilon.
$$


In code, this might be `abs(b-a) <= 2eps(m)`, which  means that the "gap" between `a` and `b` is essentially no more than $2$ floating point values from the $x$ value with the smallest $f(x)$ value.
For bracketing methods that is about as good as you can get. However, once floating point values are understood, the absolute best you can get for a bracketing interval would be


* along the way, a value `f(c)` is found which evaluates *exactly* to `0.0`

* the endpoints of the bracketing interval are *adjacent* floating point values, meaning the interval can not be bisected and `f` changes sign between the two values.


There can be problems when the stopping criteria is `abs(b-a) <= 2eps(m))` and the answer is `0.0` that require engineering around. As such, an *absolute* tolerance might be needed, one where $\lvert x_i - x_{i-1} \rvert \leq  \delta$.

For bracketing algorithms,  consideration of  $\lvert x_i - x_{i-1} \rvert$ might be all that matters, but not for algorithms like Newton's or the secant algorithm. In Newton's method the update step is $f(x_{i-1})/f'(x_{i-1})$. Naturally when $f(x_i)$ is close to $0$, the update step is small and $\lvert x_{i} - x_{i-1}\rvert = \Delta$ will be close to $0$. *However*, should $f'(x_i)$ be large, then $\Delta$ can also be small and the algorithm will possibly stop, as $x_{i} \approx x_{i-1}$---but not necessarily $x_{i} \approx \alpha$. So termination on $\Delta$ alone can be off. Checking if $f(x_{i})$ is an approximate zero---as it should be if $f$ is continuous---is  also useful to include in a stopping criteria.

However, there may never be a value with `f(x_i)` exactly `0.0`. (The value of `sin(1pi)` is not zero, for example, as `1pi` is an approximation to $\pi$, as well the `sin` of values adjacent to `float(pi)` do not produce `0.0` exactly.)


Suppose `x_i` is the closest floating point number to $\alpha$, the mathematical zero. Then the relative rounding error, $($ `x_i` $- \alpha)/\alpha$, will be a value $\delta$ with $\delta$ less than `eps()`.


How far then can `f(x_i)` be from $0 = f(\alpha)$? Consider:


$$
f(x_i) = f(x_i - \alpha + \alpha) =  f(\alpha + \alpha \cdot \delta) = f(\alpha \cdot (1 + \delta)),
$$

where $\delta = x_i/\alpha - 1$ is close to $0$ if $x_i$ converges to $\alpha$.


Assuming $f$ has a derivative, the linear approximation gives:


$$
f(x_n) \approx f(\alpha) + f'(\alpha) \cdot (\alpha\delta) = \alpha \cdot f'(\alpha) \cdot \delta
$$

So we should consider `f(x_i)` an *approximate zero* when it is on the scale of $\alpha  \cdot f'(\alpha) \cdot \delta$. That $\alpha$ factor means we consider a *relative* tolerance, $\delta$, for $f(x_i)$ based on $\lvert x_i\rvert$.

As well though, for $\alpha$ values close to $0$ this relative tolerance might be an issue, and a small absolute tolerance can be needed.


A good condition to check if `f(x_i)` is small is

* `abs(f(x_i)) <= abs(x_i) * rtol + atol`, or
* `abs(f(x_i)) <= max(abs(x_i) * rtol, atol)`


where the relative tolerance, `rtol`, would absorb an estimate for $f'(\alpha)$.



One thing to keep in mind is that the right-hand side of the rule `abs(f(x_i)) <= abs(x_i) * rtol + atol`, as a function of `x_i`, goes to `Inf` as `x_i` increases. So if `f` has `0` as an asymptote (like `e^(-x)`) for large enough `x_i`, the rule will be `true` and `x_i` could be counted as an approximate zero, despite it not being one.


A modified criteria for convergence might look like:


* stop if $\Delta$ is small and `f` is an approximate zero with some tolerances

* stop if `f` is an approximate zero with some tolerances, but be mindful that this rule can identify mathematically erroneous answers.


It is not uncommon to assign `rtol` to have a value like `sqrt(eps())` to account for accumulated floating point errors and the factor of $f'(\alpha)$, though in the `Roots` package it is set smaller by default.


### Conditioning and stability

This next part is a technical, mathematical---not practical---motivation for why we might stop when $x_i \approx x_{i-1}$ or $f(x_i) \approx 0$.

In Part III of @doi:10.1137/1.9781611977165 we find language of numerical analysis useful to  formally describe the zero-finding problem. Key concepts are errors, conditioning, and stability, which can be used to give some theoretical justification for the tolerances above.

Abstractly a *problem* is a mapping, $F$, from a domain $X$ of data to a range $Y$ of solutions. Both $X$ and $Y$ have a sense of distance given by a *norm*. A norm (denoted with $\lVert\cdot\rVert$) is a generalization of the absolute value and gives quantitative meaning to terms like small and large.

::: {.definition title="Well conditioned problem"}

A *well-conditioned* problem is one with the property that all small perturbations of $x$ lead to only small changes in $F(x)$.

:::

This sense of "small" is measured through a *condition number*.

If we let $\delta_x$ be a small perturbation of $x$ then $\delta_F = F(x + \delta_x) - F(x)$.

The *forward error* is $\lVert\delta_F\rVert = \lVert F(x+\delta_x) - F(x)\rVert$, the *relative forward error* is  $\lVert\delta_F\rVert/\lVert F\rVert = \lVert F(x+\delta_x) - F(x)\rVert/ \lVert F(x)\rVert$.

The *backward error* is $\lVert\delta_x\rVert$, the *relative backward error* is $\lVert\delta_x\rVert / \lVert x\rVert$.

 The *absolute condition number*, $\hat{\kappa}$, is the worst case of the forward error divided by the backward error, or this ratio $\lVert\delta_F\rVert/ \lVert\delta_x\rVert$, as the perturbation size shrinks to $0$.

The *relative condition number*, $\kappa$, divides $\lVert\delta_F\rVert$ by $\lVert F(x)\rVert$ and $\lVert\delta_x\rVert$ by $\lVert x\rVert$ before taking the ratio.


A *problem* is a mathematical concept, an *algorithm* the computational version. Algorithms may differ for many reasons, such as floating point errors, tolerances, etc. We use notation $\tilde{F}$ to indicate the algorithm.

The *absolute error in the algorithm* is $\lVert\tilde{F}(x) - F(x)\rVert$, the relative error divides by $\lVert F(x)\rVert$. A good algorithm would have smaller relative errors.

An algorithm is called *stable* if

$$
\frac{\lVert\tilde{F}(x) - F(\tilde{x})\rVert}{\lVert F(\tilde{x})\rVert},
$$

is *small* for *some* $\tilde{x}$ relatively near $x$, $\lVert\tilde{x}-x\rVert/\lVert x\rVert$.

::: {.definition tilte="Stable algorithm"}

A *stable* algorithm gives nearly the right answer to nearly the right question.

:::

The right answer is $F(x)$, the nearly right answer is $F(\tilde{x})$, the nearly right question is $\tilde{F}(x)$.



A related concept is an algorithm $\tilde{F}$ for a problem $F$ is *backward stable* if for each $x \in X$,

$$
\tilde{F}(x) = F(\tilde{x})
$$

for *some* $\tilde{x}$ where $\lVert\tilde{x} - x\rVert/\lVert x\rVert$ is small.

::: {.definition tilte="Backward stable algorithm"}

"A backward stable algorithm gives exactly the right answer to nearly the right question."

:::

The nearly right question is $\tilde{F}(x)$, the exactly right answer to this is $F(\tilde{x})$.


The concepts are related by Trefethen and Bao's Theorem 15.1 which says  for a backward stable algorithm the relative error $\lVert\tilde{F}(x) - F(x)\rVert/\lVert F(x)\rVert$ is small in a manner proportional to the relative condition number.

Applying this to the zero-finding we follow @doi:10.1137/1.9781611975086.

To be specific, the problem, $F$, is finding a zero of a function $f$ starting at an initial point $x_0$. The data is $(f, x_0)$, the solution is $r$ a zero of $f$.

For concreteness, take the algorithm as Newton's method. Any implementation must incorporate tolerances, so this is a computational approximation to the problem. The data is the same, but technically we use $\tilde{f}$ for the function, as any computation is dependent on machine implementations. The output is $\tilde{r}$ an *approximate* zero.

Suppose for sake of argument that $\tilde{f}(x) = f(x) + \epsilon$, $f$ has a continuous derivative, and $r$ is a root of $f$ and $\tilde{r}$ is a root of $\tilde{f}$. Then by linearization:

$$
\begin{align*}
0 &= \tilde{f}(\tilde r) \\
  &= f(r + \delta) + \epsilon\\
  &\approx f(r) + f'(r)\delta + \epsilon\\
  &= 0 + f'(r)\delta + \epsilon
\end{align*}
$$
Rearranging gives $\lVert\delta/\epsilon\rVert \approx 1/\lVert f'(r)\rVert$. But the $|\delta|/|\epsilon|$ ratio is related to the condition number:

::: {.definition title="Absolute condition number"}

The absolute condition number is $\hat{\kappa}_r = |f'(r)|^{-1}$.

:::


The error formula in Newton's method measuring the distance between the actual root and an approximation includes the derivative in the denominator, so we see large condition numbers are tied into possibly larger errors.

Now consider $g(x) = f(x) - f(\tilde{r})$. Call $f(\tilde{r})$ the residual. We have $g$ is near $f$ if the residual is small. The algorithm will solve $(g, x_0)$ with $\tilde{r}$, so with a small residual an exact solution to an approximate question will be found. Driscoll and Braun state

::: {.relationship title="Backward error and residual"}

The backward error in a root estimate is equal to the residual.

:::

Practically these two observations lead to

* If there is a large condition number, it may not be possible to find an approximate root near the real root.

* A tolerance in an algorithm should consider both the size of $x_{n} - x_{n-1}$ and the residual $f(x_n)$.

For the first observation, the example of Wilkinson's polynomial is often used where $f(x) = (x-1)\cdot(x-2)\cdot \cdots\cdot(x-20)$. When expanded this function has exactness issues of typical floating point values, the condition number is large and some of the roots found are quite different from the mathematical values.


The second observation follows from $f(x_n)$ monitoring the backward error and the product of the condition number and the backward error monitoring the forward error. This product is on the order of $|f(x_n)/f'(x_n)|$ or $|x_{n+1} - x_n|$.


## Questions


###### Question


Let `f(x) = tanh(x)` (the hyperbolic tangent) and `fp(x) = sech(x)^2`, its derivative.


Does *Newton's* method (using `Roots.Newton()`) converge starting at `1.0`?


```{julia}
#| hold: true
#| echo: false
yesnoq("yes")
```

Does *Newton's* method (using `Roots.Newton()`) converge starting at `1.3`?


```{julia}
#| hold: true
#| echo: false
yesnoq("no")
```

Does the secant method (using `Roots.Secant()`) converge starting at `1.3`? (a second starting value will automatically be chosen, if not directly passed in.)


```{julia}
#| hold: true
#| echo: false
yesnoq("yes")
```

###### Question


For the function `f(x) = x^5 - x - 1` both Newton's method and the secant method will converge to the one root when started from `1.0`. Using `verbose=true` as an argument to `find_zero`, (e.g., `find_zero(f, x0, Roots.Secant(), verbose=true)`) how many *more* steps does the secant method need to converge?


```{julia}
#| hold: true
#| echo: false
numericq(2)
```

Do the two methods converge to the exact same value?


```{julia}
#| hold: true
#| echo: false
yesnoq("yes")
```

###### Question


Let `f(x) = exp(x) - x^4` and `x0=8.0`. How many steps (iterations) does it take for the secant method to converge using the default tolerances?


```{julia}
#| hold: true
#| echo: false
numericq(10, 1)
```

###### Question


Let `f(x) = exp(x) - x^4` and a starting bracket be `x0 = [8, 9]`. Then calling `find_zero(f,x0, verbose=true)` will show that 48 steps are needed for exact bisection to converge. What about with the `Roots.Brent()` algorithm, which uses inverse quadratic steps when it can?


It takes how many steps?


```{julia}
#| hold: true
#| echo: false
numericq(36, 1)
```

The `Roots.A42()` method uses inverse cubic interpolation, as possible, how many steps does this method take to converge?


```{julia}
#| hold: true
#| echo: false
numericq(7, 1)
```

The large difference is due to how the tolerances are set within `Roots`. The Brent method gets pretty close in a few steps, but takes a much longer time to get close enough for the default tolerances.


###### Question


Consider this crazy function defined by:


```{julia}
#| eval: false
f(x) = cos(100*x)-4*erf(30*x-10)
```

(The `erf` function is the [error function](https://en.wikipedia.org/wiki/Error_function) and is in the `SpecialFunctions` package loaded with `CalculusWithJulia`.)


Make a plot over the interval $[-3,3]$ to see why it is called "crazy".


Does `find_zero` find a zero to this function starting from $0$?


```{julia}
#| hold: true
#| echo: false
yesnoq("yes")
```

If so, what is the value?


```{julia}
#| hold: true
#| echo: false
f(x) = cos(100*x)-4*erf(30*x-10)
val = find_zero(f, 0)
numericq(val)
```

If not, what is the reason?


```{julia}
#| hold: true
#| echo: false
choices = [
"The zero is a simple zero",
"The zero is not a simple zero",
"The function oscillates too much to rely on the tangent line approximation far from the zero",
"We can find an answer"
]
answer = 4
radioq(choices, answer, keep_order=true)
```

Does `find_zero` find a zero to this function starting from $0.175$?


```{julia}
#| hold: true
#| echo: false
yesnoq(false)
```

If so, what is the value?


```{julia}
#| hold: true
#| echo: false
numericq(-999.999)
```

If not, what is the reason?


```{julia}
#| hold: true
#| echo: false
choices = [
"The zero is a simple zero",
"The zero is not a simple zero",
"The function oscillates too much to rely on the tangent line approximations far from the zero",
"We can find an answer"
]
answer = 3
radioq(choices, answer, keep_order=true)
```
