Newton's method and bisection in Python

If you are in Year 5 of the applied sciences liceo (the Italian academic high school with a science and computing focus), in the first term computer science reaches numerical methods: finding the zeros of a function with the computer, when the equation cannot be solved with a formula. The written test or lab assignment on this block usually falls six to nine weeks from now, between late November and December. There is time, and the topic is one of the best in the programme, because for the first time the code you write gives a result you can check with a calculator.

Every year in my lessons I see the same two obstacles: understanding why the method converges, and writing a loop with a sensible stopping criterion. I always use the same example, the simplest there is: $f(x)=x^2-2$, whose positive root is $\sqrt{2}\approx 1.41421356$. We already know the answer, so we can see whether the method finds it.

If Python is still a bit new to you, start from programming in Python from scratch. To think about loops and conditions before the code, you need flowcharts and pseudocode. For another example of an algorithm with a cost comparison, see sorting algorithms.

The bisection method

The idea comes from the intermediate value theorem: if $f$ is continuous on $[a,b]$ and $f(a)$ and $f(b)$ have opposite signs, there is at least one root in the interval. You take the midpoint $m=(a+b)/2$ and keep the half where the sign changes. Repeat.

For $f(x)=x^2-2$ on $[1,2]$: $f(1)=-1$ and $f(2)=2$, opposite signs, so we can start.

After 5 steps the root lies in $[1.40625;1.4375]$: the maximum error is $0.03125$. At every step the interval halves, so after $n$ steps the width is $(b-a)/2^n$. To get an error below $10^{-6}$ from an interval of width 1 you need $2^n>10^6$, that is $n=20$ (indeed $2^{20}=1,048,576$).

The Newton–Raphson method

Newton’s method uses the derivative: from the point $x_n$ you draw the tangent to the graph and take the point where it meets the $x$-axis. The formula is:

For $f(x)=x^2-2$ the derivative is $f^{\prime }(x)=2x$, and the formula becomes $x_{n+1}=x_n-\frac{x_n^2-2}{2x_n}$. Let us start from $x_0=1.5$:

• $x_1=1.5-\frac{0.25}{3}=1.41666\dots$
• $x_2=1.41666\dots -\frac{0.006944\dots }{2.8333\dots }=1.41421568\dots$
• $x_3=1.41421356237\dots$

After only 3 steps there are 11 correct digits. A useful observation for the test: the number of correct digits roughly doubles at every step (error: about $2.5\cdot 10^{-3}$, then $2\cdot 10^{-6}$, then $10^{-12}$). Bisection, by contrast, gains just one bit (a factor of 2) per step.

The Python code

I tested these functions with python3 on the example above.

def f(x):
    return x * x - 2

def df(x):
    return 2 * x

def bisection(f, a, b, tol=1e-6, max_iter=100):
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    for i in range(1, max_iter + 1):
        m = (a + b) / 2
        if f(m) == 0 or (b - a) / 2 < tol:
            return m, i
        if f(a) * f(m) < 0:
            b = m
        else:
            a = m
    return (a + b) / 2, max_iter

def newton(f, df, x0, tol=1e-10, max_iter=50):
    x = x0
    for i in range(1, max_iter + 1):
        d = df(x)
        if d == 0:
            raise ZeroDivisionError("zero derivative")
        x_new = x - f(x) / d
        if abs(x_new - x) < tol:
            return x_new, i
        x = x_new
    return x, max_iter

print(bisection(f, 1, 2))
print(newton(f, df, 1.5))

Output: bisection gives about 1.4142141 in 20 iterations, Newton 1.4142135623730951 in 4 iterations. With a tolerance of $10^{-10}$ bisection needs 34. For the same precision, Newton does less than 15% of the work.

A few choices you should be able to explain in the lab:

• the functions receive f as a parameter, so the same code works for any equation;
• max_iter avoids an infinite loop if the method does not converge;
• the initial sign check in bisection protects against intervals that do not contain a root;
• in newton the d == 0 check avoids division by zero.

Stopping criteria

A numerical loop must know when to stop. There are three criteria, often combined:

1. Width of the interval (bisection): $(b-a)/2<tol$.
2. Difference between approximations: $|x_{n+1}-x_n|<tol$.
3. Residual: $|f(x_n)|<tol$. Careful: if the function is very flat near the root, a small residual does not guarantee an $x$ close to the root.

Always add a maximum number of iterations: it is the last safety net.

Comparison and limits

In practice: Newton is faster, bisection is more robust. A typical mistake is to start Newton from $x_0=0$ for $x^2-2$. The derivative $2x$ is 0 and the code raises the exception. With $x_0=0.1$ the first step sends the point to $x_1=0.1-\frac{-1.99}{0.2}=10.05$, and more steps are needed to get back near the root: the method works, but the starting point matters.

If you want a level beyond high school, where error and convergence are also studied formally, you will find it in numerical methods for the exam. For the calculus behind Newton’s method, it is worth revising derivatives explained well.

A complete exercise

Find $\sqrt{3}$ as the positive root of $f(x)=x^2-3$, using Newton from $x_0=2$. The derivative is $2x$ and the step is $x_{n+1}=x_n-\frac{x_n^2-3}{2x_n}$.

• $x_1=2-\frac{1}{4}=1.75$
• $x_2=1.75-\frac{0.0625}{3.5}=1.7321428\dots$
• $x_3=1.7320508100\dots$

The exact value is $\sqrt{3}=1.7320508075\dots$: at the third step the error is already below $10^{-8}$. How many iterations would bisection have used on $[1,2]$ for the same tolerance? You need $2^n>10^8$, that is $n=27$ ($2^{27}=134,217,728$). The comparison 3 against 27 is the answer many teachers want to hear about the speed of convergence.

A warning about the code: if you are asked to modify it for $x^2-3$, only f and df change, not the two method functions. If you had to touch bisection or newton, something is not general.

Test yourself

Before a test, redo by hand three iterations of Newton for $f(x)=x^2-3$ from $x_0=2$ and then check with the code. Then try the free practice test on high school tests. If you want to go over your lab code together, we can do it in the free 30-minute introductory call with me.