A note on Steffensen's method for solving equations

Steffensen's method is based on Newton's iteration for solving a non-linear equation \(f(x)=0\):

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

Newton's method can fail to work in a number of ways, but when it does work it displays qudratic convergence; the number of correct signifcant figures roughly doubling at each step. However, it also has the disadvntage of needing to compute the derivative as well as the function. This may be difficult for some functions.

Steffensen's idea was to use the quotient approximation of the derivative:

\[ f'(x)\approx\frac{f(x+h)-f(x)}{h} \]

when \(h\) is small, and since we trying to solve \(f(x)=0\), we may assume that in the neighbourhood of the solution \(f(x)\) is itself small, so can be used for \(h\). This means we can write

\[ f'(x)\approx\frac{f(x+f(x))-f(x)}{f(x)} \]

which leads to the following version of Newton's method:

\[ x\leftarrow x-\frac{f(x)^2}{f(x+f(x))-f(x)}. \]

This is a neat idea, and in fact when it works it converges almost as fast as Newton's method. However, it is very senstive to the starting value. For example, suppose we want to find the value of \(W(10)\), where \(W(x)\) is Lambert's \(W\) function; the inverse of \(y=xe^x\). Finding \(W(10)\) then means solving the equation

\[ xe^x-10=0. \]

Newton's method uses the iteration

\[ x\leftarrow x-\frac{xe^x-10}{e^x(x+1)} \]

and with a positive starting value not too big will converge; the first 50 places of the solution are:

\[ 1.74552800274069938307430126487538991153528812908093 \]

Staring with \(x=2\) will produce over 1000 correct decimal places in 12 steps.

If we apply Steffensen's method starting with \(x=2\) we'll see values that wobble about 1.9 for ages 1.9 before converging to the wrong value. Newton's method will work for almost any value (although the larger the initial value, the long the iterations take to "settle down"); Steffensen's method will only work when the initial value is close to 1.7.

A slight improvement

Using the "central" approximation of the derivative:

\[ f'(x)\approx\frac{f(x+h)-f(x-h)}{2h} \]

makes a considerable difference; this leads to the iteration

\[ x\leftarrow x-\frac{2f(x)^2}{f(x+f(x))-f(x-f(x))} \]

This does however require the computation of three function values, rather than just the original two. A slightly faster version of the above is

\[ x\leftarrow x-\frac{f(x)^2}{f(x+\frac{1}{2}f(x))-f(x-\frac{1}{2}f(x))}. \]