The Newton method for finding roots

In science and engineering we constantly need to solve equations like $f(x) = 0$. A point $x^*$ with $f(x^*) = 0$ is called a root (or zero) of $f$.

The Newton method is a classic iterative algorithm that finds such a root by following the direction in which the function is pointing: at the current guess it builds the tangent line and follows it down to the $x$-axis to obtain the next, hopefully better, guess.

Everything below the sliders runs in your browser as WebAssembly — no Julia server, no install.

begin
    using PlutoUI, WasmMakie
end

The idea

Suppose we have a guess $x_n$ for the root. The tangent line at $x_n$ has slope $f'(x_n)$, so it crosses zero at

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$$

Repeating this gives a sequence $x_0, x_1, x_2, \dots$ that converges very quickly (quadratically) once it is close to a root.

The original lecture used automatic differentiation (ForwardDiff.jl) to get $f'$. WebAssembly can't run that machinery, so here we instead supply hand-written analytic derivatives for a small menu of example functions — the mathematics is exactly the same.

Pick an example function

choice	$f(x)$	$f'(x)$	a root
1	$x^2 - 2$	$2x$	$\sqrt 2 \approx 1.41421$
2	$x^3 - x$	$3x^2 - 1$	$0,\ \pm 1$
3	$\cos x - x$	$-\sin x - 1$	$\approx 0.73909$

example function = 1

starting point x₀ = 2.0

number of iterations n = 5

Visualising the steps

The curve is $f$. The horizontal line is the $x$-axis (where the root lives). The vertical ticks mark each iterate $x_0, x_1, \dots$ as Newton's method homes in on the root.

let
    # sample the curve with an integer loop (StepRangeLen iteration is not
    # wasm-compilable yet) over a fixed window
    xs = Float64[]
    ys = Float64[]
    lo = -3.0
    hi = 3.0
    steps = 240
    for k in 0:steps
        t = lo + (hi - lo) * k / steps
        push!(xs, t)
        push!(ys, f(t, which))
    end

    # the Newton iterates (x positions converging on the root)
    seq = newton_sequence(which, x0, n)

    # proven-minimal WasmMakie API (bare Axis + bare lines!, like the known-good
    # figure island): the curve f(x), the x-axis line, and a vertical tick at
    # each Newton iterate so you can watch them home in on the root.
    fig = Figure(size = (480, 320))
    ax = Axis(fig[1, 1])
    lines!(ax, xs, ys)                  # the curve f(x)
    lines!(ax, [lo, hi], [0.0, 0.0])    # the x-axis (where the root lives)
    for v in seq
        lines!(ax, [v, v], [-0.5, 0.5]) # a tick marking this iterate
    end
    fig
end

Final estimate of the root: 1.4142135623730951

Residual |f(x)| at that estimate: 4.440892098500626e-16

(this shrinks toward 0 as you raise the iteration count n)

The example function f(x) selected by the slider (which = 1, 2 or 3).

"The example function `f(x)` selected by the slider (`which` = 1, 2 or 3)."
function f(x::Float64, which::Int64)
    if which == 1
        return x*x - 2.0
    elseif which == 2
        return x*x*x - x
    else
        return cos(x) - x
    end
end

fprime

The hand-written analytic derivative fprime(x) matching f.

"The hand-written analytic derivative `fprime(x)` matching `f`."
function fprime(x::Float64, which::Int64)
    if which == 1
        return 2.0*x
    elseif which == 2
        return 3.0*x*x - 1.0
    else
        return -sin(x) - 1.0
    end
end

The iteration

newton_sequence runs the Newton update $x_{n+1} = x_n - f(x_n)/f'(x_n)$ a fixed number of times and records every iterate, so we can watch it converge.

newton_sequence

Return the Newton iterates x₀, x₁, …, xₙ for example which starting at x0.

"Return the Newton iterates x₀, x₁, …, xₙ for example `which` starting at `x0`."
function newton_sequence(which::Int64, x0::Float64, n::Int64)
    xs = Vector{Float64}(undef, n + 1)
    x = x0
    xs[1] = x
    for i in 1:n
        d = fprime(x, which)
        # guard against a flat tangent (division by zero)
        if d == 0.0
            xs[i + 1] = x
        else
            x = x - f(x, which) / d
            xs[i + 1] = x
        end
    end
    return xs
end

"2.0"

"1.5"

"1.4166666666666667"

"1.4142156862745099"

"1.4142135623746899"

"1.4142135623730951"

# shown as strings so the Float64 list renders identically in-browser
# (whole-number iterates like 2.0 print as "2.0", matching Julia)
iterates = string.(newton_sequence(which, x0, n))

newton_root

The final Newton estimate of the root after n iterations.

"The final Newton estimate of the root after `n` iterations."
function newton_root(which::Int64, x0::Float64, n::Int64)
    xs = newton_sequence(which, x0, n)
    return xs[end]
end

1.4142135623730951

root_estimate = newton_root(which, x0, n)

How well did it converge?

residual is $|f(x_n)|$ at the final iterate — it should shrink toward 0 as you increase the number of iterations n (try the slider!).

residual

The absolute residual |f(x)| at the final Newton iterate.

"The absolute residual |f(x)| at the final Newton iterate."
function residual(which::Int64, x0::Float64, n::Int64)
    x = newton_root(which, x0, n)
    return abs(f(x, which))
end

4.440892098500626e-16

res = residual(which, x0, n)

Appendix

Newton's method converges quadratically near a simple root — the number of correct digits roughly doubles each step — but it can fail if it lands on a point where $f'(x) = 0$ (a flat tangent), or wander off if it starts far from any root. Try example 2 ($x^3 - x$) with different starting points to see both fast convergence and the method jumping between the three roots $0, \pm 1$.