32  Other zero-finding algorithms

This section uses these add-on packages:

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.

32.1 Other methods

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

32.1.1 Estimating the derivative

Sidi 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\).

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

pprime_n2 (generic function with 1 method)

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.

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₃

4-element Vector{Float64}:
  1.0
 -0.7853981633974482
  0.31440769948782343
 -0.043893642219757556

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:

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

1.0298665293222589

32.1.2 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\).

@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

I (generic function with 1 method)

We use \(k=2\) and see what comes:

I(2)

\(\frac{- x₁ y₂ + x₂ y₁}{y₁ - y₂}\)

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

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)

\(0\)

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.

u = lambdify(I(3), (x[1:3]..., y[1:3]...))

Callable function with variables (:x₁, :x₂, :x₃, :y₁, :y₂, :y₃)

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

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

3-element Vector{Float64}:
  1.0
 -0.7853981633974482
  0.31440769948782343

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

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

(1.0695976825969624, -0.054321738106927975)

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

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

1.0298665293222589

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

An inverse cubic interpolation is utilized by Alefeld, Potra, and Shi 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.

32.1.3 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 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.

32.1.4 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:

@syms xᵢ xᵢ₋₁ yᵢ yᵢ₋₁ x
only(solve(yᵢ / (x - xᵢ) ~ -yᵢ₋₁ / (xᵢ₋₁ - x), x))

\(\frac{- xᵢ yᵢ₋₁ + xᵢ₋₁ yᵢ}{yᵢ - yᵢ₋₁}\)

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.

Figure 32.1 show that some function shapes result in one end point being fixed as the algorithm progresses which can lead to linear convergence.

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

Modified regula falsi

Figure 32.2 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.

Figure 32.2: 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 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 Figure 32.3.

Figure 32.3: 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 Figure 32.1.

NoteHybrid 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.

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

(0, 1)

The Steffensen method needs a nearby estimate:

find_zero(f, x0, Roots.Steffensen())   # 5 iterations, 12 function evaluations

0.8745112203760901

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:

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

0.8745112203760902

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

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

0.8745112203760902

Finally, we compare regula falsi variants:

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

0.8745112203760902

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.

32.2 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.

32.2.1 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 Trefethen and Bau (2022) 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.

DefinitionWell 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\).

DefinitionDefinition

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.

DefinitionDefinition

“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 Driscoll and Braun (2017).

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:

DefinitionAbsolute 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

RelationshipBackward 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|\).

32.3 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?

Select an item

Yes

No

Does Newton’s method (using Roots.Newton()) converge starting at 1.3?

Select an item

Yes

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.)

Select an item

Yes

No

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?

Do the two methods converge to the exact same value?

Select an item

Yes

No

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?

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?

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

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:

f(x) = cos(100*x)-4*erf(30*x-10)

(The erf function is the 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\)?

Select an item

Yes

No

If so, what is the value?

If not, what is the reason?

Select an item

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

Does find_zero find a zero to this function starting from \(0.175\)?

Select an item

Yes

No

If so, what is the value?

If not, what is the reason?

Select an item

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

Driscoll, Tobin A., and Richard J. Braun. 2017. Fundamentals of Numerical Computation. Julia adaptation. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9781611975086.

Trefethen, Lloyd N., and David Bau. 2022. Numerical Linear Algebra, Twenty-Fifth Anniversary Edition. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9781611977165.