Newton’s method with 10 lines of Python

danielhomola Blog 6 Comments

I’m starting a new series of blog posts, called “XY in less than 10 lines of Python“. This first one is about Newton’s method, which is an old numerical approximation technique that could be used to find the roots of complex polynomials and any differentiable function.

Let’s say we have a complicated polynomial:

f(x)=6x^5-5x^4-4x^3+3x^2

and we want to find its roots. Unfortunately we know from the Galois theory that there is no formula for a 5th degree polynomial so we’ll have to use numeric methods. Wait a second. What’s going on? We all learned the quadratic formula in school, and there are formulas for cubic and quartic polynomials, but Galois proved that no such “root-finding” formula exist for fifth or higher degree polynomials, that uses only the usual algebraic operations (addition, subtraction, multiplication, division) and application of radicals (square roots, cube roots, etc).

Some nice guys pointed out on reddit that I didn’t quite get the theory right. Sorry about that, I’m no mathematician by any means. It turns out that this polynomial could be factored into x^2(x-1)(6x^2 + x - 3) and solved with traditional cubic formula. Also the theorem I referred to is the Abel-Ruffini Theorem and it only applies to the solution to the general polynomial of degree five or greater. Nonetheless the example is still valid, and demonstrates how would you apply Newton’s method, to any polynomial, so let’s crack on.

So in these cases we have to resort to numeric linear approximation. A simple way would be to use the intermediate value theorem, which states that if f(x) is continuous on [a,b] and f(a) < y< f(b), then there is an x between a and b so that f(x)=y.

We could exploit this by looking for an x_1 where f(x_1)>0 and an x_2 where f(x_2)<0, and then we could be certain that f(x)=0 must be between x_1 and x_2. So we could  check x_3=\frac{x_2-x_1}{2}, and find that f(x_3) will be positive, then continue withx_4=\frac{x_2-x_3}{2} which would be negative and our proposed x_n would become closer and closer to f(x_n)=0. This method however can be pretty slow, so Newton devised a better way to speed things up (when it works).

If you look at the figure below, you’ll see the plot of our polynomial. It has three roots at 0, 1, and somewhere in-between. So how do we find these?

In Newton’s method we take a random point f(x_0), then draw a tangent line through x_0, f(x_0), using the derivative f'(x_0). The point x_1 where this tangent line crosses the x axis will become the next proposal we check. We calculate the tangent line at f'(x_1) and find x_2. We carry on, and as we do |f(x_n)| \to 0, or in other words we can make our approximation as close to zero as we want, provided we are willing to continue with the number crunching. Using the fact that the slope of tangent (by the definition of derivatives) at x_n is f'(x_n) we can derive the formula for x_{n+1}, i.e. where the tangent crosses the x axis:

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

Newton's method

With all this in mind it’s easy to write an algorithm that approximates $f(x)=0$ with arbitrary error \varepsilon. Obviously in this example depending on where we start x_0 we might find different roots.

There you go, we’are done in 10 lines (9 without the blank line), even less without the print statements. So in order to use this, we need two functions, f(x) and f'(x). For this polynomial these are:

Now we can simply find the three roots with:

All of the above code, and some additional comparison test with the scipy.optimize.newton method can be found in this Gist. And don’t forget, if you find it too much trouble differentiating your functions, just use SymPy, I wrote about it here.

Newton’s method is pretty powerful but there could be problems with the speed of convergence, and awfully wrong initial guesses might make it not even converge ever, see here. Nonetheless I hope you found this relatively useful.. Let me know in the comments.

danielhomolaNewton’s method with 10 lines of Python
  • Analytics2Go

    Love it!

  • Very cool. Congrats!

  • Ericracia

    Hey, How do you put code in your blog’s posts like this one? Is it a plug-in?

  • Силантий

    Why function dx returns the absolute value of the difference between zero and the value of the function? If you can just return the module value of the function. It’s the same thing.