Newton method for interval root finding

Download the notebook for this tutorial.


This page will explain you how to generalize the traditional Newton method for root-finding to deal with intervals. Basically, here you will get a taste of what happens under the hood in the IntervalRootFinding when the Newton method is used. First, let us import the packages we will need

using IntervalArithmetic, ForwardDiff, Plots

Review: traditional Newton method

It is well known from high-school or undergraduate studies that given an initial guess x0x_0 for the root of a function ff, the root can be found iteratively with the upgrade rule

xk+1=xkf(xk)f(xk), x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)},

which is known as Newton method. This has also the following geoemtrical interpretation:

  1. Draw the tangent to ff passing through the point (xk,f(xk))(x_k, f(x_k))

  2. The new point xk+1x_{k+1} will be the root of the tangent line.

Let us demonstrate this with a simple example suppose we want to find the root of the function f(x)=atan(x2)f(x) = atan(\frac{x}{2}) using x0=2.5x_0=2.5 as initial guess. The following animation shows the Newton method in action. Note that the derivative can be computed using the ForwardDiff package

f(x) = atan(x/2)
x0 = 2.5
y0 = f(x0)
d(f, x) = ForwardDiff.derivative(f, x)

xk = 2.5
yk = f(xk)

anim = @animate for i in 1:5
    xnew = xk - yk/d(f, xk)
    ynew = f(xnew)
    scatter([xk], [0]; markersize=5, markershape=:circle, c=:black, label="\$ x_k\$")
    plot!([xk, xk], [0, yk], linestyle=:dash, color=:black, label=nothing)
    plot!([xnew, xk], [0, yk], color=:black, label=nothing)
    scatter!([xnew], [0], color=:blue, markersize=5, label="\$x_{k+1}\$")
    plot!(f, -10, 10, label="\$f(x)\$", c=:red, legend=:right, legendfont=12, xlims=(-5, 5))
    global xk, yk = xnew, ynew

@show xk, yk
(xk, yk) = (-3.911818029136288e-8, -1.9559090145681437e-8)

After 55 iterations the absolute error is already in the magnitude 10810^{-8}. If the initial guess x0x_0 is close enough to the correct root, then the Newton method can give a fast and accurate solution. However, it has no information about the error. This motives to find a generalized interval version.

Interval Newton method in 1D

Suppose we want to find the roots of a function ff over the interval XX. Let us take as initial guess the the midpoint of XX m(X)m(X). The core idea is now to consider all possible slopes a tangent line to ff in XX can have, instead of just the tangent at m(X)m(X). This means we define the Newton operator

Nf(X)=m(X)f(m(X))f(X) N_f(X) = m(X) - \frac{f(m(X))}{f'(X)}

where f(X)f'(X) takes into acount all possible slopes the tangent can have in XX. As an example, consider the function f(x)=atan(x2)f(x)=atan(\frac{x}{2}) with starting interval X=[1,5]X=[-1,5]. Observe the following figure

X = -1..5
N(f, X) = mid(X) - f(mid(X))/d(f, X)
X = -1..5

Nx = N(f, X)
Xnext = Nx ∩ X

x0 = mid(X)
y0 = f(x0)
offset = 5e-2
plot([X.lo, X.hi], [-offset, -offset], label="\$X\$", c=:green, lw=2)
scatter!([x0], [-offset]; markersize=5, markershape=:circle, c=:black, label=nothing)
plot!([x0, x0], [-offset, y0], linestyle=:dash, color=:black, label=nothing)

deriv = d(f, X)
for m in range(deriv.lo, deriv.hi, length=100) # plot cone
    plot!([x0-y0/m, x0], [offset, y0], color=:gray, label=nothing, alpha=0.5)
plot!([Nx.lo, Nx.hi], [offset, offset], label="\$N(X)\$", c=:blue, lw=2)
plot!([Xnext.lo, Xnext.hi], [0, 0], label="\$N(x) \\cap X\$", c=:black, lw=3)
plot!(f, -10, 10, label="\$f(x)\$", c=:red, legend=:right)
plot!(legend=:bottomright, legendfont=12)

The green interval represents our starting interval. The gray area is the cone containing all lines passing through (m(X), f(m(X))) with slope in f(X)f'(X). The blue interval is Nf(X)N_f(X). Now, since we new the root has to been contained both in XX and in Nf(X)N_f(X), we can conclude that the root must lie in XNf(X)X \cap N_f(X), which is highlighted in dark in the figure. As you can notice, we squeezed the starting interval to a smaller one, which is guaranteed to contain the root. We are now ready to present the iteration of the interval newton method

Xk+1=(m(Xk)m(Xk)f(Xk))Xk X_{k+1} = (m(X_k) - \frac{m(X_k)}{f'(X_k)}) \cap X_k

Now the following important theorem holds

Theorem if 0f(X)0 \notin f'(X), then the intervals XkX_k form a nested sequence squeezing quadratically to the root of ff.

Note that the condition 0f(X)0 \notin f'(X) implies that either ff has a unique root in XX or no roots at all. An equivalent expression for this theorem is the following:

  • If Nf(X)XN_f(X)\subseteq X, then XX has a unique root

  • if Nf(X)X=N_f(X)\cap X=\emptyset, then ff has no root in XX.

The following animation shows the Newton method in action with our example from before

function newton_iteration(f, df, X)
    m = @interval mid(X)
    return (m - f(m)/df(X)) ∩ X

xk = X
df(x) = d(f, x)
anim = @animate for i in 1:4
    zerobox = IntervalBox(xk, f(X))
    plot!(f, -1, 1, c=:red, xlims=(-1,1), legend=nothing)
    global xk = newton_iteration(f, df, xk)
    @show xk
xk = [-1, 0.4292036733]
xk = [-0.001913860913, 0.06895721471]
xk = [-3.670747462e-05, 3.138504583e-06]
xk = [-3.940431152e-16, 5.259983123e-15]

As you can see, after only 44 iterations the interval has already shrunk to an interval of width 101510^{-15}, which is guaranteed to contain the true value of the root.

We have seen that the Newton method is very powerful when 0f(X)0\notin f'(X), but what if this condition is not met? If 0f(X)0\in f'(X), then division by zero occurs. This means that using traditional interval arithmetic the ratio f(m(X))f(X)\frac{f(m(X))}{f'(X)} will be [,][-\infty, \infty], which is not very useful. To overcome this problem we exploit extended division, suppose 0f(Xk)=Y=[Y,Y]0\in f'(X_k)=Y=[\underline{Y}, \overline{Y}] , then we consider the two intervals Y1=[Y,0]Y_1=[\underline{Y},0] and Y2=[0,Y]Y_2=[0,\overline{Y}]. Then we apply the update formula using both intervals obtaining the two intervals

Xk+1a=(m(Xk)m(Xk)Y1)Xk X_{k+1}^a=\left(m(X_k) - \frac{m(X_k)}{Y_1}\right) \cap X_k Xk+1b=(m(Xk)m(Xk)Y2)Xk X_{k+1}^b=\left(m(X_k) - \frac{m(X_k)}{Y_2}\right) \cap X_k

and then repeat Newton iteration with extended division to both interals. Hence, in the worst-case scenario we will double the number of intervals to track. If at some point we have Xk=X_k=\emptyset, then we can conclude it has no roots and it can bo thrown away. If Nf(Xk)XkN_f(X_k)\subseteq X_k, then the root is unique and we do not need to split the interval with extended division anymore.

Generalization to higher dimensions

The Interval Newton method can be generalized to higher dimensions replacing the derivative with the Jacobian. Suppose you want to solve the system F(x)=0F(\mathbf{x})=\mathbf{0}, where F:RnRnF:\mathbb{R}^n\rightarrow\mathbb{R}^n. Now the update formula for interval Newton is

Xk+1=m(Xk)J(Xk)1F(m(Xk)), X_{k+1} = m(X_k) - J(X_k)^{-1}F(m(X_k)),

where XkX_k is now an interval box.