Interpolation-based methods - Fundamentals of Numerical Computation

From a practical standpoint, one of the biggest drawbacks of Newton’s method is the requirement to supply $f'$ in Function 4.3.2. 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

Julia

MATLAB

Python

Graphical interpretation of the secant method

We return to finding a root of the equation $x e^x=2$ .

f(x) = x * exp(x) - 2;

plot(f, 0.25, 1.25, label="function",
    xlabel=L"x", ylabel=L"y", legend=:topleft)

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.

x₁ = 1;
y₁ = f(x₁);
x₂ = 0.5;
y₂ = f(x₂);
scatter!([x₁, x₂], [y₁, y₂], label="initial points",
    title="Two initial values")

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.

m₂ = (y₂ - y₁) / (x₂ - x₁)
secant = x -> y₂ + m₂ * (x - x₂)
plot!(secant, 0.25, 1.25, label="secant line", l=:dash, color=:black,
    title="Secant line")

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

x₃ = x₂ - y₂ / m₂
@show y₃ = f(x₃)
scatter!([x₃], [0], label="root of secant", title="First iteration")

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.

m₃ = (y₃ - y₂) / (x₃ - x₂)
x₄ = x₃ - y₃ / m₃

Graphical interpretation of the secant method

We return to finding a root of the equation $x e^x=2$ .

f = @(x) x .* exp(x) - 2;
clf, fplot(f, [0.25, 1.25])
set(gca, 'ygrid', 'on')  
xlabel('x'), ylabel('y')    
title('Objective function')

x1 = 1;    y1 = f(x1);
x2 = 0.5;  y2 = f(x2);
hold on, scatter([x1, x2], [y1, y2])
title('Two initial values')

slope2 = (y2 - y1) / (x2 - x1);
secant2 = @(x) y2 + slope2 * (x - x2);

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

fplot(secant2,[0.25, 1.25],'k--')
x3 = x2 - y2 / slope2;
y3 = f(x3)
scatter(x3, 0)
title('Next value')

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.

slope2 = (y3 - y2) / (x3 - x2);
x4 = x3 - y3 / slope2;
y4 = f(x4)

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

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

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

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

The example in Demo 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 iteration formula.

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

Algorithm 4.4.2 (secant)

Julia

MATLAB

Python

Secant method

secant.m

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
function x = secant(f,x1,x2)
% SECANT   Secant method for a scalar equation.
% Input:
%   f        objective function 
%   x1,x2    initial root approximations
% Output       
%   x        vector of root approximations (last is best)

% Operating parameters.
funtol = 100*eps;  xtol = 100*eps;  maxiter = 40;

x = [x1 x2];
dx = Inf;  y1 = f(x1);
k = 2;  y2 = f(x2);

while (abs(dx) > xtol) && (abs(y2) > funtol) && (k < maxiter)
    dx = -y2 * (x(k)-x(k-1)) / (y2-y1);   % secant step
    x(k+1) = x(k) + dx;
    
    k = k+1;
    y1 = y2;    % current f-value becomes the old one next time
    y2 = f(x(k));
end

if k==maxiter
    warning('Maximum number of iterations reached.')
end

Secant method

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]

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 tedious 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

Julia

MATLAB

Python

Convergence of the secant method

We check the convergence of the secant method from Demo 4.4.1. Again we will use extended precision to get a longer sequence than double precision allows.

f(x) = x * exp(x) - 2
x = FNC.secant(f, BigFloat(1), BigFloat(0.5), xtol=1e-80, ftol=1e-80);

We don’t know the exact root, so we use the last value as a proxy.

r = x[end]

Here is the sequence of errors.

ϵ = @. Float64(r - x[1:end-1])

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.

[log(abs(ϵ[k+1])) / log(abs(ϵ[k])) for k in 1:length(ϵ)-1]

As expected, this settles in at around 1.618.

Convergence of the secant method

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

f = @(x) x .* exp(x) - 2;
x = secant(f, 1, 0.5);

We don’t know the exact root, so we use fzero to get a good proxy.

r = fzero(f, 1);

Here is the sequence of errors.

format short e
err = r - x(1:end-1)'

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));
ratios = logerr(2:end) ./ logerr(1:end-1)

As expected, this settles in at around 1.618.

Convergence of the secant method

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

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

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)

Here is the sequence of errors.

err = r - x
print(err)

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

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

Julia

MATLAB

Python

Inverse quadratic interpolation

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

f(x) = x + cos(10 * x)
interval = [0.5, 1.5]

plot(f, interval..., label="Function", legend=:bottomright,
    grid=:y, ylim=[-0.1, 3], xlabel=L"x", ylabel=L"y")

We choose three values to get the iteration started.

x = [0.8, 1.2, 1]
y = @. f(x)
scatter!(x, y, label="initial points")

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 = Polynomials.fit(x, y, 2)      # interpolating polynomial
plot!(x -> q(x), interval..., l=:dash, label="interpolant")

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

:::{card}
By giving two functions in the plot call, we get the parametric plot $(q(y),y)$ as a function of $y$.

plot(f, interval..., label="Function",
    legend=:bottomright, grid=:y, xlabel=L"x", ylabel=L"y")
scatter!(x, y, label="initial points")

q = Polynomials.fit(y, x, 2)       # interpolating polynomial
plot!(y -> q(y), y -> y, -0.1, 2.6, l=:dash, label="inverse interpolant")

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

q(0)

Let’s restart the process with BigFloat numbers to get a convergent sequence.

x = BigFloat.([8, 12, 10]) / 10
y = @. f(x)

for k = 3:12
    q = Polynomials.fit(y[k-2:k], x[k-2:k], 2)
    push!(x, q(0))
    push!(y, f(x[k+1]))
end

println("residual = $(f(x[end]))")

As far as our current precision is concerned, we have an exact root.

r = x[end]
logerr = @. log(Float64(abs(r - x[1:end-1])))
[logerr[k+1] / logerr[k] for k in 1:length(logerr)-1]

The convergence is probably superlinear at a rate of $\alpha=1.8$ or greater.

Inverse quadratic interpolation

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

f = @(x) x + cos(10 * x);
interval = [0.5, 1.5];
clf, fplot(f, interval)
set(gca, 'ygrid', 'on'), axis(axis)   
title('Objective function')    
xlabel('x'), ylabel('y')    
r = fzero(f, 1)

We choose three values to get the iteration started.

x = [0.8, 1.2, 1]';
y = f(x);
hold on, scatter(x, y)
title('Three initial points')

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.

c = polyfit(x, y, 2);    % coefficients of interpolant
q = @(x) polyval(c, x);
fplot(q, interval, '--')
title('Parabola model')

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

:::{card}
By giving two functions in the `fplot` call, we get the parametric plot $(q(y),y)$ as a function of $y$.

cla, fplot(f, interval)
scatter(x, y)     
c = polyfit(y, x, 2);    % coefficients of interpolating polynomial
q = @(y) polyval(c, y);
fplot(q, @(y) y, ylim,'--')    % plot x=q(y), y=y
title('Sideways parabola')

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

x = [x; q(0)];
y = [y; f(x(end))]

We repeat the process a few more times.

for k = 4:8
    c = polyfit(y(k-2:k), x(k-2:k), 2);
    x(k+1) = polyval(c, 0);
    y(k+1) = f(x(k+1));
end
disp('final residual:')
y(end)

Here is the sequence of errors.

format short e
err = x - r

The convergence is probably superlinear:

logerr = log(abs(err));
ratios = logerr(2:end) ./ logerr(1:end-1)

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)

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}")

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)

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

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) ⌨ Use nlsolve with ftol=1e-15 to find a reference value for each root.

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

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

$x^2=e^{-x}$ , over $[-2,2]$
$2x = \tan x$ , over $[-0.2,1.4]$
$e^{x+1}=2+x$ , over $[-2,2]$
⌨ Use a plot to approximately locate all the roots of $f(x)=x^{-2}-\sin(x)$ in the interval $[0.5,10]$ . Then find a pair of initial points for each root such that Function 4.4.2 converges to that root.
✍ Show analytically that the secant method converges in one step for a linear function, regardless of the initialization.
✍ In general, the secant method formula (4.4.2) cannot be applied if $x_{k}=x_{k-1}$ . However, suppose that $f(x)=ax^2+bx+c$ for constants $a$ , $b$ , and $c$ . Show that in this case the formula can be simplified to one that is well defined when $x_{k}=x_{k-1}$ . Then show that the resulting $x_{k+1}$ is the same as the result of one step of Newton’s method applied to $f$ at $x_k$ .
✍ Let $f(x)=x^2$ . Show that if $(1/x_1)$ and $(1/x_2)$ are positive integers, and the secant iteration is applied, then the sequence $1/x_1,1/x_2,1/x_3,\ldots$ is a Fibonacci sequence, i.e., satisfying $x_{k+1}=x_k+x_{k-1}$ .
✍ Provide the details that show how to derive (4.4.3) from (4.4.2).
⌨ Write a function iqi(f,x₁,x₂,x₃) that performs inverse quadratic interpolation for finding a root of $f$ , given three initial estimates. To find the quadratic polynomial $q(y)$ passing through the three most recent points, use fit. Test your function on the function in Exercise 1 from this section.

Preface

Newton’s method

Preface

Newton for nonlinear systems