Householder's methods

These are a class of root-finding methods; that is, for the numerical solution of a single nonlinear equation, developed by Alston Scott Householder in 1970. They may be considered a generalisation of the well known Newton-Raphson method (also known more simply as Newton’s method) defined by

\[ x\leftarrow x-\frac{f(x)}{f'(x)}. \]

where the equation to be solved is \(f(x)=0\).

From a starting value \(x_0\) a sequence of iterates can be generated by

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

As is well known, Newton’s method exhibits quadratic convergence; that is, if the sequence of iterates converges to a root value \(r\), then the limit

\[ \lim_{n\to\infty}\frac{x_{n+1}-r}{(x_n-r)^2} \]

is finite. This means, in effect, that the number of correct decimal places doubles at each step. Householder’s method for a rate of convergence \(d+1\) is defined by

\[ x\leftarrow x-d\frac{(1/f)^{(d-1)}(x)}{(1/f)^{(d)}(x)}. \]

We show how this definition can be rewritten in terms of ratios of derivatives, by using Python and its symbolic toolbox SymPy.

We start by defining some variables and functions.

from sympy import *
x = Symbol('x')
f = Function('f')(x)

Now we can define the first Householder formula, with \(d=1\):

d = 1
H1 = x + d*diff(1/f,x,d-1)/diff(1/f,x,d)

\[ x-\frac{f(x)}{\frac{d}{dx}f(x)} \]

which is Newton’s formula. Now for \(d=2\):

d = 2
H2 = x + d*diff(1/f,x,d-1)/diff(1/f,x,d)

\[ x - \frac{2 \frac{d}{d x} f{\left (x \right )}}{- \frac{d^{2}}{d x^{2}} f{\left (x \right )} + \frac{2 \left(\frac{d}{d x} f{\left (x \right )}\right)^{2}}{f{\left (x \right )}}} \]

This is a mighty messy formula, but it can be greatly simplified by using ratios of derivatives defined by

\[ r_k=\frac{f^{(d-1}(x)}{f^{(d)}(x)} \] This means that \[ r_1=\frac{f}{f'},\quad r_2=\frac{f'}{f^{\prime\prime}} \] To make the subsitution into the current expression above, we can use the substitutions \[ f^{\prime\prime}=f'/r_2,\quad f'=f/r_1 \] to be done sequentially (first defining the new symbols)

r_1,r_2,r_3 = symbols('r_1,r_2,r_3')
H2r = H2s.subs([(Derivative(f,x,2), Derivative(f,x)/r_2), (Derivative(f,x), f/r_1)]).simplify()

\[ -\frac{2r_1r_1}{r_1-2r_2} \] Dividing the top and bottom by \(2r_2\) produces the formulation \[ \frac{r_1}{1-\displaystyle{\frac{r_1}{2r_2}}} \] and so Householder’s method for \(d=2\) is defined by the recurrence \[ x\leftarrow x-\frac{r_1}{1-\displaystyle{\frac{r_1}{2r_2}}}. \] This is known as Halley’s method, after Edmond Halley, also known for his comet. This method has been called the most often rediscovered iteration formula in the literature.

It would exhibit cubic convergence, which means that the number of correct figures roughly triples at each step.

Apply the same sequence of steps for \(d=3\), and including the subsitution \[ f^{\prime\prime\prime} = f^{\prime\prime}/r_3 \] produces the fourth order formula \[ x\leftarrow x-\frac{3 r_{1} r_{3} \left(2r_{2} - r_{1}\right)}{r_{1}^{2} - 6 r_{1} r_{3} + 6 r_{2} r_{3}} \]

comments powered by Disqus