Higher order Newton-like methods for root finding

There are variations on Newton's root finding method that use higher derivatives and converge faster. Alston Householder developed a sequence of such methods of the form

where the superscripts in parentheses indicate derivatives. When n = 2, Householder's method reduces to Newton's method.

Once Newton's method is close enough to a root, the error is squared at each iteration. Said another way, the number of correct decimal places roughly doubles. With the Householder method of order n, the number of correct decimal places roughly increases by a factor of n.

The Householder method H_n can be written as a rational function of f and its derivatives up to order n-1, though the formulas are complicated.

Maybe it's time for Householder's methods to be more widely used: with automatic differentiation [1], the complexity of computing higher order derivatives is not the drawback it was when these computations were done solely by hand.

Demonstration

We'll illustrate the convergence of the Householder methods by finding roots of the quintic polynomial

f(x) = x⁵ + ax + b.

I used Mathematica to calculate expressions for the Householder methods and export the results as code. Mathematica does not have a function to export Python code, but it does have a function CForm to export expressions as C code. I started to use this, converting exponents from C to Python form with search and replace. Then I tried using FortranForm instead, and that was much easier because Fortran and Python both represent powers with a ** operator.

def f(x, a, b): return x**5 + a*x + bdef H2(x, a, b): return x - (b + a*x + x**5)/(a + 5*x**4) def H3(x, a, b): num = (a + 5*x**4)*(b + a*x + x**5) den = a**2 - 10*b*x**3 + 15*x**8 return x - num/dendef H4(x, a, b): num = (b + a*x + x**5)*(-a**2 + 10*b*x**3 - 15*x**8) den = a**3 + 10*b**2*x**2 + 5*a**2*x**4 - 80*b*x**7 - 25*a*x**8 + 35*x**12 return x + num/den

To set values of a and b, I arbitrarily chose a = 2 and solved for b so that I is a root. That way we have a known value of the root for comparison. I used x = 4 as my starting guess at the root.

I first used -log₁₀|x - I| to measure the number of correct decimals. That didn't quite work because the method would converge exactly (to the limits of machine precision) which led to taking the log of zero. So I modified my accuracy function to the following.

def decimals(x): d = abs(x - pi) if d != 0: return( round(-log10(d),2) ) else: return "inf"

Then I tested each of the Householder methods with the following.

epsilon = 1e-14for H in [H2, H3, H4]: x = 4 for i in range(1, 10): if abs(f(x,a,b)) > epsilon: x = H(x,a,b) print(i, x, decimals(x))

Here are the results in tabular form. I replaced all the decimals after the first incorrect decimal with underscores to make it easier to see the convergence.

|---+-------------------+--------|| i | x | digits ||---+-------------------+--------|| 1 | 3.4______________ | 0.53 || 2 | 3.18_____________ | 1.33 || 3 | 3.142____________ | 2.87 || 4 | 3.141593_________ | 5.93 || 5 | 3.141592653590___ | 12.07 || 6 | 3.141592653589793 | inf ||---+-------------------+--------|| 1 | 3.2______________ | 1.11 || 2 | 3.1416___________ | 4.03 || 3 | 3.1415926535899__ | 12.79 || 4 | 3.141592653589793 | inf ||---+-------------------+--------|| 1 | 3.15_____________ | 1.84 || 2 | 3.141592654______ | 8.85 || 3 | 3.141592653589793 | inf ||---+-------------------+--------|

Note that, as promised, the number of correct decimal places for H_n increases by roughly a factor of n with each iteration.

Update: See this post where the symbolic calculations are done in SymPy rather than Mathematica. Everything stays in Python, and there's no need to edit the code for the Householder methods.

Related posts

• Number of real roots in an interval
• Bessel function crossings

[1] Note that in this post I used Mathematica to do symbolic differentiation, not automatic differentiation. Automatic differentiation differs from both symbolic and numerical differentiation. See an explanation here.