from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

4.4.Interpolation-based methods¶

From a practical standpoint, one of the biggest drawbacks of Newton’s method is the requirement to supply $f'$ in Function 4.3.1. It is both a programming inconvenience and a step that requires computational time. We can avoid using $f'$ , however, by making a simple but easily overlooked observation:

Let’s call this the principle of approximate approximation.

In the Newton context, the principle of approximate approximation begins with the observation that the use of $f'$ is linked to the construction of a linear approximation $q(x)$ equivalent to a tangent line. The root of $q(x)$ is used to define the next iterate in the sequence. We can avoid calculating the value of $f'$ by choosing a different linear approximation.

Example 4.4.1 (Graphical interpretation of the secant method)

f = lambda x: x * exp(x) - 2
xx = linspace(0.25, 1.25, 400)

fig, ax = subplots()
ax.plot(xx, f(xx), label="function")
ax.set_xlabel("$x$")
ax.set_ylabel("$f(x)$")
ax.grid();

From the graph, it’s clear that there is a root near $x=1$ . To be more precise, there is a root in the interval $[0.5,1]$ . So let us take the endpoints of that interval as two initial approximations.

x1 = 1
y1 = f(x1)
x2 = 0.5
y2 = f(x2)
ax.plot([x1, x2], [y1, y2], "ko", label="initial points")
ax.legend()
fig

Instead of constructing the tangent line by evaluating the derivative, we can construct a linear model function by drawing the line between the two points $\bigl(x_1,f(x_1)\bigr)$ and $\bigl(x_2,f(x_2)\bigr)$ . This is called a secant line.

slope2 = (y2 - y1) / (x2 - x1)
secant2 = lambda x: y2 + slope2 * (x - x2)
ax.plot(xx, secant2(xx), "--", label="secant line")
ax.legend()
fig

As before, the next root estimate in the iteration is the root of this linear model.

x3 = x2 - y2 / slope2
ax.plot(x3, 0, "o", label="root of secant")
y3 = f(x3)
print(y3)
ax.legend()
fig

-0.17768144843679456

For the next linear model, we use the line through the two most recent points. The next iterate is the root of that secant line, and so on.

slope3 = (y3 - y2) / (x3 - x2)
x4 = x3 - y3 / slope3
print(f(x4))

0.05718067370113333

The example in Example 4.4.1 demonstrates the secant method. In the secant method, one finds the root of the linear approximation through the two most recent root estimates. That is, given previous approximations $x_1,\ldots,x_k$ , define the linear model function as the line through $\bigl(x_{k-1},f(x_{k-1})\bigr)$ and $\bigl(x_k,f(x_k)\bigr)$ :

q(x) = f(x_k) + \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}(x-x_k).

(4.4.1)

Solving $q(x_{k+1})=0$ for $x_{k+1}$ gives the following iteration formula.

Our implementation of the secant method is given in Function 4.4.1.

Algorithm 4.4.1 (secant)

secant.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def secant(f, x1, x2):
    """
    secant(f, x1, x2)

    Use the secant method to find a root of f starting from x1 and x2. Returns a
    vector of root estimates.
    """
    # Operating parameters.
    eps = np.finfo(float).eps
    funtol = 100 * eps
    xtol = 100 * eps
    maxiter = 40

    x = np.zeros(maxiter)
    x[:2] = [x1, x2]
    y1 = f(x1)
    y2 = 100
    dx = np.inf  # for initial pass below
    k = 1

    while (abs(dx) > xtol) and (abs(y2) > funtol) and (k < maxiter):
        y2 = f(x[k])
        dx = -y2 * (x[k] - x[k - 1]) / (y2 - y1)  # secant step
        x[k + 1] = x[k] + dx  # new estimate

        k = k + 1
        y1 = y2  # current f-value becomes the old one next time

    if k == maxiter:
        warnings.warn("Maximum number of iterations reached.")
    return x[:k+1]

About the code

Because we want to observe the convergence of the method, Function 4.4.1 stores and returns the entire sequence of root estimates. However, only the most recent two are needed by the iterative formula. This is demonstrated by the use of y₁ and y₂ for the two most recent values of $f$ .

4.4.1Convergence¶

Graphically, a secant line usually looks like a less accurate model of $f$ than the tangent line. How will that affect the convergence?

As before, let $\epsilon_k = x_k-r$ be the errors in the successive root approximations, and assume that $r$ is a simple root. If the initial errors are small, then a lengthy but straightforward Taylor expansion shows that, to lowest order,

\epsilon_{k+1} \approx \frac{1}{2}\frac{f''(r)}{f'(r)} \epsilon_k \epsilon_{k-1}.

(4.4.3)

If we make an educated guess that

\epsilon_{k+1} = c (\epsilon_k)^\alpha, \quad \epsilon_k = c (\epsilon_{k-1})^\alpha, \ldots, \qquad \alpha>0,

(4.4.4)

then (4.4.3) becomes

\left[ \epsilon_{k-1}^{\alpha} \right]^{\,\alpha} \approx C \epsilon_{k-1}^{\alpha+1}

(4.4.5)

for an unknown constant $C$ . Treating the approximation as an equality, this becomes solvable if and only if the exponents match, i.e., $\alpha^2 = \alpha+1$ . The only positive root of this equation is the golden ratio,

\alpha = \frac{1+\sqrt{5}}{2} \approx 1.618.

(4.4.6)

Hence, the errors in the secant method converge like $\epsilon_{k+1} = c (\epsilon_k)^\alpha$ for $1<\alpha<2$ .

Quadratic convergence is a particular case of superlinear convergence. Roughly speaking, we expect

\begin{align*} \log |\epsilon_{k+1}| & \approx \alpha (\log |\epsilon_k|) + \log L, \\ \frac{\log |\epsilon_{k+1}|}{\log |\epsilon_k|} & \approx \alpha + \frac{\log L}{\log |\epsilon_k|} \to \alpha, \end{align*}

(4.4.8)

as $k\to\infty$ .

Example 4.4.2 (Convergence of the secant method)

We check the convergence of the secant method from Example 4.4.1.

f = lambda x: x * exp(x) - 2
x = FNC.secant(f, 1, 0.5)
print(x)

[1.         0.5        0.81037177 0.86563193 0.85217802 0.85260123
 0.8526055  0.8526055  0.8526055 ]

We don’t know the exact root, so we use root_scalar to get a substitute.

from scipy.optimize import root_scalar
r = root_scalar(f, bracket=[0.5, 1]).root
print(r)

0.8526055020137254

Here is the sequence of errors.

err = r - x
print(err)

[-1.47394498e-01  3.52605502e-01  4.22337271e-02 -1.30264253e-02
  4.27479941e-04  4.26991559e-06 -1.40547707e-09  4.55191440e-15
  0.00000000e+00]

It’s not easy to see the convergence rate by staring at these numbers. We can use (4.4.8) to try to expose the superlinear convergence rate.

logerr = log(abs(err))
for i in range(len(err) - 2):
    print(logerr[i+1] / logerr[i])

0.5444386280277934
3.0358017547194565
1.3716940021941457
1.7871469297607958
1.593780475055435
1.6485786717829443
1.62014464365765

/var/folders/0m/fy_f4_rx5xv68sdkrpcm26r00000gq/T/ipykernel_67680/1382054591.py:1: RuntimeWarning: divide by zero encountered in log
  logerr = log(abs(err))

As expected, this settles in at around 1.618.

In terms of the error as a function of the iteration number $k$ , the secant method converges at a rate strictly between linear and quadratic, which is slower than Newton’s method. But error versus iteration count may not be the best means of comparison.

Often we analyze rootfinding methods by assuming that the bulk of computing time is spent evaluating the user-defined functions $f$ and $f'$ . (Our simple examples and exercises mostly don’t support this assumption, but many practical applications do.) In this light, we see that Newton’s method requires two evaluations, $f(x_k)$ and $f'(x_k)$ , for each iteration. The secant method, on the other hand, while it uses the two function values $f(x_k)$ and $f(x_{k-1})$ at each iteration, only needs to compute a single new one. Note that Function 4.4.1 keeps track of one previous function value rather than recomputing it.

Now suppose that $|\epsilon_k|=\delta$ . Roughly speaking, two units of work (i.e., function evaluations) in Newton’s method brings us to an error of $\delta^2$ . If one spreads out the improvement in the error evenly across the two steps, using

\delta^2 = \bigl( \delta^{\sqrt{2}} \bigr)^{\!\sqrt{2}},

(4.4.9)

it seems reasonable to say that the rate of convergence in Newton per function evaluation is $\sqrt{2}\approx 1.41$ . This is actually less than the comparable rate of about 1.62 for the secant method.

Not only is the secant method easier to apply than Newton’s method in practice, it’s also more efficient—a rare win-win!

4.4.2Inverse interpolation¶

At each iteration, the secant method constructs a linear model function that interpolates the two most recently found points on the graph of $f$ . Two points determine a straight line, so this seems like a sensible choice. But as the iteration progresses, why use only the two most recent points? What would it mean to use more of them?

If we interpolate through three points by a polynomial, we get a unique quadratic function. Unfortunately, a parabola may have zero, one, or two crossings of the $x$ -axis, potentially leaving some doubt as to how to define the next root estimate. On the other hand, if we turn a parabola on its side, we get a graph that intersects the $x$ -axis exactly once, which is ideal for defining the next root estimate.

This leads to the idea of defining $q(y)$ as the quadratic interpolant to the points $(y_{k-2},x_{k-2})$ , $(y_{k-1},x_{k-1})$ , and $(y_k,x_k)$ , where $y_i=f(x_i)$ for all $i$ , and setting $x_{k+1}=q(0)$ . The process defined in this way (given three initial estimates) is called inverse quadratic interpolation. Rather than deriving lengthy formulas for it here, we demonstrate how to perform inverse quadratic interpolation using fit to perform the interpolation step.

Example 4.4.3 (Inverse quadratic interpolation)

Here we look for a root of $x+\cos(10x)$ that is close to 1.

f = lambda x: x + cos(10 * x)
xx = linspace(0.5, 1.5, 400)
fig, ax = subplots()
ax.plot(xx, f(xx), label="function")
ax.grid()
xlabel("$x$"), ylabel("$y$")
fig

We choose three values to get the iteration started.

x = array([0.8, 1.2, 1])
y = f(x)
ax.plot(x, y, "ko", label="initial points")
ax.legend()
fig

If we were using forward interpolation, we would ask for the polynomial interpolant of $y$ as a function of $x$ . But that parabola has no real roots.

q = poly1d(polyfit(x, y, 2))  # interpolating polynomial
ax.plot(xx, q(xx), "--", label="interpolant")
ax.set_ylim(-0.1, 3), ax.legend()
fig

To do inverse interpolation, we swap the roles of $x$ and $y$ in the interpolation.

plot(xx, f(xx), label="function")
plot(x, y, "ko", label="initial points")

q = poly1d(polyfit(y, x, 2))  # inverse interpolating polynomial
yy = linspace(-0.1, 2.6, 400)
plot(q(yy), yy, "--", label="inverse interpolant")

grid(), xlabel("$x$"), ylabel("$y$")
legend();

We seek the value of $x$ that makes $y$ zero. This means evaluating $q$ at zero.

x = hstack([x, q(0)])
y = hstack([y, f(x[-1])])
print("x:", x, "\ny:", y)

x: [0.8        1.2        1.         1.10398139] 
y: [0.65449997 2.04385396 0.16092847 1.14820652]

We repeat the process a few more times.

for k in range(6):
    q = poly1d(polyfit(y[-3:], x[-3:], 2))
    x = hstack([x, q(0)])
    y = hstack([y, f(x[-1])])
print(f"final residual is {y[-1]:.2e}")

final residual is 1.47e-14

Here is the sequence of errors.

from scipy.optimize import root_scalar
r = root_scalar(f, bracket=[0.9, 1]).root
err = x - r
print(err)

[-1.67888402e-01  2.32111598e-01  3.21115982e-02  1.36092984e-01
  1.53473435e-02  3.26831473e-03  4.61743614e-04  6.29584770e-06
  3.43903706e-09  6.37268016e-14]

The error seems to be superlinear, but subquadratic:

logerr = log(abs(err))
for i in range(len(err) - 1):
    print(logerr[i+1] / logerr[i])

0.818477543986624
2.354297090262193
0.5800188694773643
2.0942526256501486
1.3702985897217017
1.341928284536839
1.5592238835587917
1.6273123181772837
1.5591161372818163

4.4.3Bracketing¶

Like Newton’s method, the secant and inverse quadratic interpolation methods cannot guarantee convergence. One final new idea is needed to make a (nearly) foolproof algorithm.

If $f$ is continuous on the interval $[a,b]$ and $f(a)f(b)<0$ —that is, $f$ changes sign on the interval—then $f$ must have at least one root in the interval, due to the Intermediate Value Theorem from calculus. If we come up with a new root estimate $c\in(a,b)$ , then whatever sign $f(c)$ is, it is different from the sign at one of the endpoints. (Of course, if $f(c)$ is zero, we are done!) So either $[a,c]$ or $[c,b]$ is guaranteed to have a root too, and in this way we can maintain not just individual estimates but an interval that always contains a root.

The best algorithms blend the use of fast-converging methods with the guarantee provided by a bracket. For example, say that an iteration begins with a bracketing interval. Make a list of the inverse quadratic estimate, the secant estimate, and the midpoint of the current interval, and pick the first member of the list that lies within the current interval. Replace the interval with the bracketing subinterval, and start a new iteration. This is the idea behind Brent’s method, which is a very successful rootfinding algorithm.

4.4.4Exercises¶

For each of Exercises 1–3, do the following steps.

(a) ✍ Rewrite the equation into the standard form for rootfinding, $f(x) = 0$ .

(b) ⌨ Make a plot of $f$ over the given interval and determine how many roots lie in the interval.

(c) ⌨ Find a reference value for each root.

(d) ⌨ Determine a bracketing interval for each root. Then use Function 4.4.1, starting with the endpoints of the bracketing interval, to find each root.

(e) ⌨ For one of the roots, use the errors in the resulting sequence to determine numerically whether the convergence is apparently between linear and quadratic.