Steffensen's method

In

, whereas the more familiar Newton's method also converges quadratically, but requires derivatives and the secant method does not require derivatives but also converges less quickly than quadratically.

Steffensen's method has the drawback that it requires two function evaluations per step, whereas the secant method requires only one evaluation per step, so it is not necessarily most efficient in terms of

Steffensen's method can be derived as an adaptation of Aitken's delta-squared process applied to fixed-point iteration. Viewed in this way, Steffensen's method naturally generalizes to efficient fixed-point calculation in general Banach spaces, whenever fixed points are guaranteed to exist and fixed-point iteration is guaranteed to converge, although possibly slowly, by the Banach fixed-point theorem.

The simplest form of the formula for Steffensen's method occurs when it is used to find a

${\displaystyle f}$; that is, to find the real value ${\displaystyle \ x_{\star }\ }$ that satisfies ${\displaystyle \ f(x_{\star })=0~.}$ Near the solution ${\displaystyle \ x_{\star }\ ,}$ the derivative of the function, ${\displaystyle \ f'\ ,}$ needs to either exactly or very nearly satisfy ${\displaystyle -1<f'(x_{\star })<0~.}$^[b] For some functions, Steffensen's method can work even if this condition is not met, but in such a case, the starting value ${\displaystyle \ x_{0}\ }$ must be very close to the actual solution ${\displaystyle \ x_{\star }\ ,}$ and convergence to the solution may be slow. Adjustment of the size of the method's intermediate step, mentioned later, can improve convergence in some of these cases.

Given an adequate starting value ${\displaystyle \ x_{0}\ ,}$ a sequence of values ${\displaystyle \ x_{0},\ \ x_{1},\ x_{2},\ \dots ,\ x_{n},\ \dots \ }$ can be generated using the formula below. When it works, each value in the sequence is much closer to the solution ${\displaystyle \ x_{\star }\ }$ than the prior value. The value ${\displaystyle \ x_{n}\ }$ from the current step generates the value ${\displaystyle \ x_{n+1}\ }$ for the next step, via the formula^[1]

${\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{g(x_{n})}}}$

for ${\displaystyle \ n=0,1,2,3,...\ ,}$ where the slope function ${\displaystyle \ g(x)\ }$ is a composite of the original function ${\displaystyle \ f\ }$ given by the formula

${\displaystyle g(x)={\frac {f{\bigl (}x+f(x){\bigr )}}{f(x)}}-1}$

or perhaps more clearly,

${\displaystyle g(x)={\frac {f(x+h)-f(x)}{h}}\qquad \approx \quad {\frac {\operatorname {d} f(x)}{\operatorname {d} x}}\equiv f'(x),}$

where ${\displaystyle \ h\equiv f(x)\ }$ is a step-size between the last iteration point, ${\displaystyle \ x\ ,}$ and an auxiliary point located at ${\displaystyle \ x+h~.}$

Technically, the function ${\displaystyle \ g\ }$ is called the first-order

of ${\displaystyle \ f\ }$ between those two points^[c] Practically, it is the averaged value of the slope ${\displaystyle f'}$ of the function ${\displaystyle \ f\ }$ between the last sequence point ${\displaystyle \left(x,y\right)={\bigl (}x_{n},f\left(x_{n}\right){\bigr )}}$ and the auxiliary point at ${\displaystyle \ {\bigl (}x,y{\bigr )}={\bigl (}x_{n}+h,f\left(x_{n}+h\right){\bigr )}\ ,}$ with the size of the intermediate step (and its direction) given by ${\displaystyle \ h=f(x_{n})~.}$

Because the value of ${\displaystyle \ g\ }$ is an approximation for ${\displaystyle \ f'\ ,}$ its value can optionally be checked to see if it meets the condition ${\displaystyle \ -1<g<0\ ,}$ which is required to guarantee convergence of Steffensen's algorithm. Although slight non-conformance may not necessarily be dire, any large departure from the condition warns that Steffensen's method is liable to fail, and temporary use of some fallback algorithm is warranted (e.g. the more robust

).

It is only for the purpose of finding ${\displaystyle \ h\ }$ for this auxiliary point that the value of the function ${\displaystyle \ f\ }$ must fulfill the requirement that ${\displaystyle \ -1<f'(x_{\star })<0~.}$[b] For all other parts of the calculation, Steffensen's method only requires the function ${\displaystyle \ f\ }$ to be continuous and to actually have a nearby solution.^[1] Several modest modifications of the step ${\displaystyle \ h\ }$ used in the formula for the slope ${\displaystyle \ g\ }$ exist, such as multiplying it by ⁠ 1 /2⁠ or ⁠ 3 /4⁠, to accommodate functions ${\displaystyle \ f\ }$ that do not quite meet the requirement.

The main advantage of Steffensen's method is that it has

– that is, both methods find roots to an equation ${\displaystyle \ f\ }$ just as "quickly". In this case, quickly means that for both methods, the number of correct digits in the answer doubles with each step. But the formula for Newton's method requires evaluation of the function's derivative ${\displaystyle \ f'\ }$ as well as the function ${\displaystyle \ f\ ,}$ while Steffensen's method only requires ${\displaystyle \ f\ }$ itself. This is important when the derivative is not easily or efficiently available.

The price for the quick convergence is the double function evaluation: Both ${\displaystyle \ f(x_{n})\ }$ and ${\displaystyle \ f(x_{n}+h)\ }$ must be calculated, which might be time-consuming if ${\displaystyle \ f\ }$ is complicated. For comparison, both regula falsi and the secant method only need one function evaluation per step. The secant method increases the number of correct digits by "only" a factor of roughly   1.6 per step,  but one can do twice as many steps of the secant method within a given time. Since the secant method can carry out twice as many steps in the same time as Steffensen's method,^[d] in practical use the secant method actually converges faster than Steffensen's method, when both algorithms succeed: The secant method achieves a factor of about   (1.6)²   ≈   2.6 times   as many digits for every two steps (two function evaluations), compared to Steffensen's factor of   2   for every one step (two function evaluations).

Similar to most other iterative root-finding algorithms, the crucial weakness in Steffensen's method is choosing a "sufficiently close" starting value ${\displaystyle \ x_{0}~.}$ If the value of ${\displaystyle \ x_{0}\ }$ is not "close enough" to the actual solution ${\displaystyle \ x_{\star }\ ,}$ the method may fail, and the sequence of values ${\displaystyle \ x_{0},\,x_{1},\,x_{2},\,x_{3},\,\dots \ }$ may either erratically flip-flop between two (or more) extremes, or diverge to infinity, or both.

The version of Steffensen's method implemented in the MATLAB code shown below can be found using Aitken's delta-squared process for convergence acceleration. To compare the following formulae to the formulae in the section above, notice that ${\displaystyle x_{n}=p-p_{n}}$. This method assumes starting with a linearly convergent sequence and increases the rate of convergence of that sequence. If the signs of ${\displaystyle p_{n},\,p_{n+1},\,p_{n+2}}$ agree and ${\displaystyle p_{n}}$ is "sufficiently close" to the desired limit of the sequence ${\displaystyle p}$, then we can assume

${\displaystyle {\frac {p_{n+1}-p}{p_{n}-p}}\approx {\frac {p_{n+2}-p}{p_{n+1}-p}},}$

so that

${\displaystyle (p_{n+2}-2p_{n+1}+p_{n})p\approx p_{n+2}p_{n}-p_{n+1}^{2}.}$

Solving for the desired limit of the sequence ${\displaystyle p}$ gives:

${\displaystyle p\approx {\frac {p_{n+2}p_{n}-p_{n+1}^{2}}{p_{n+2}-2p_{n+1}+p_{n}}}}$

${\displaystyle =~{\frac {\,(\,p_{n}^{2}+p_{n}\,p_{n+2}-2\,p_{n}\,p_{n+1}\,)-(\,p_{n}^{2}-2\,p_{n}\,p_{n+1}+p_{n+1}^{2}\,)\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}}$

${\displaystyle =p_{n}-{\frac {(p_{n+1}-p_{n})^{2}}{p_{n+2}-2p_{n+1}+p_{n}}},}$

which results in the more rapidly convergent sequence:

${\displaystyle p\approx p_{n+3}=p_{n}-{\frac {(p_{n+1}-p_{n})^{2}}{p_{n+2}-2p_{n+1}+p_{n}}}.}$

Here is the source for an implementation of Steffensen's Method in MATLAB.

function Steffensen(f, p0, tol)
% This function takes as inputs: a fixed point iteration function, f, 
% and initial guess to the fixed point, p0, and a tolerance, tol.
% The fixed point iteration function is assumed to be input as an
% inline function. 
% This function will calculate and return the fixed point, p, 
% that makes the expression f(x) = p true to within the desired 
% tolerance, tol.

format compact   % This shortens the output.
format long      % This prints more decimal places.

for i = 1:1000   % get ready to do a large, but finite, number of iterations.
                 % This is so that if the method fails to converge, we won't
                 % be stuck in an infinite loop.
    p1 = f(p0) + p0;  % calculate the next two guesses for the fixed point.
    p2 = f(p1) + p1;
    p = p0-(p1-p0)^2/(p2-2*p1+p0) % use Aitken's delta squared method to
                                  % find a better approximation to p0.
    if abs(p - p0) < tol  % test to see if we are within tolerance.
        break             % if we are, stop the iterations, we have our answer.
    end
    p0 = p;               % update p0 for the next iteration.
end

if abs(p - p0) > tol      % If we fail to meet the tolerance, we output a
                          % message of failure.
    'failed to converge in 1000 iterations.'
end

Here is the source for an implementation of Steffensen's method in Python.

from typing import Callable, Iterator
Func = Callable[[float], float, float]

def g(f: Func, x: float, fx: float) -> Func:
    """First-order divided difference function.

    Arguments:
        f: Function input to g
        x: Point at which to evaluate g
        fx: Function f evaluated at x 
    """
    return lambda x: f(x + fx) / fx - 1

def steff(f: Func, x: float, tol: float) -> Iterator[float]:
    """Steffenson algorithm for finding roots.

    This recursive generator yields the x_{n+1} value first then, when the generator iterates,
    it yields x_{n+2} from the next level of recursion.

    Arguments:
        f: Function whose root we are searching for
        x: Starting value upon first call, each level n that the function recurses x is x_n
    """
    n = 0

    while True:

        if n > 1000:
           print( "failed to converge in 1000 itterations" )
           break
        else:
           n = n + 1

        fx = f(x)

        if abs(fx) < tol:
            break
        else:
            gx = g(f, x, fx)(x)
            x = x - fx / gx    # Update to x_{n+1}
            yield x            # Yield value

Steffensen's method can also be used to find an input ${\displaystyle \ x=x_{\star }\ }$ for a different kind of function ${\displaystyle \ F\ }$ that produces output the same as its input: ${\displaystyle \ x_{\star }=F(x_{\star })\ }$ for the special value ${\displaystyle \ x_{\star }~.}$ Solutions like ${\displaystyle \ x_{\star }\ }$ are called

.

Momentarily ignoring the issues of a more general

real numbers for the sake of an example: To re-orient the reader to the earlier section, a simple toy model

fixed-point function, ${\displaystyle \ {\tilde {F}}\ ,}$ using any root function ${\displaystyle \ f\ ,}$ can be made with ${\displaystyle \ {\tilde {F}}(x)=x+\varepsilon \ f(x)~.}$ Here ${\displaystyle \ \varepsilon \ }$ is a constant with the appropriate sign that is small enough in magnitude to make ${\displaystyle \ {\tilde {F}}\ }$ stable under iteration, but large enough for the

non-linearity

of the function ${\displaystyle \ f\ }$ to be appreciable.

This method for finding fixed points of a real-valued function has been generalized for functions ${\displaystyle \ F:X\to X\ }$ that map a Banach space ${\displaystyle \ X\ }$ onto itself or even more generally ${\displaystyle \ F:X\to Y\ }$ that map from one Banach space ${\displaystyle X}$ into another Banach space ${\displaystyle \ Y~.}$ The generalized method assumes that a

linear operators

${\displaystyle \ {\bigl \{}\ G(u,v):u,v\in X\ {\bigr \}}\ }$ associated with ${\displaystyle \ u\ }$ and ${\displaystyle \ v\ }$ can be devised that (locally) satisfies the condition^[2]

${\displaystyle F\left(u\right)-F\left(v\right)=G\left(u,v\right)\ {\bigl (}\ u-v\ {\bigr )}\quad }$

1

The operator ${\displaystyle \ G\ }$ is roughly equivalent to a

vector arguments

${\displaystyle \ u\ }$ and ${\displaystyle \ v~.}$ Refer again back to the simple function ${\displaystyle \ f\ ,}$ given in the first section, where the function merely takes in and puts out real numbers: There, the function ${\displaystyle \ g\ }$ is a

divided difference

. In the generalized form here, the operator ${\displaystyle \ G\ }$ is the analogue of a divided difference for use in the Banach space.

If division is possible in the Banach space, then the linear operator ${\displaystyle \ G\ }$ can be obtained from

${\displaystyle G\left(u,v\right)={\bigl [}\ F\left(u\right)-F\left(v\right)\ {\bigr ]}\ {\bigl (}\ u-v\ {\bigr )}^{-1}\ ,}$

which may provide some insight: Expressed in this way, the linear operator ${\displaystyle \ G\ }$ can be more easily seen to be an elaborate version of the

divided difference

${\displaystyle \ g\ }$ discussed in the first section, above. The quotient form is shown here for orientation only; it is not required per se. Note also that division within the Banach space is not necessary for the elaborated Steffensen's method to be viable; the only requirement is that the operator ${\displaystyle \ G\ }$ satisfy (1).

Steffensen's method is then very similar to the Newton's method, except that it uses the divided difference ${\displaystyle \ G{\bigl (}F\left(x\right),x{\bigr )}\ }$ instead of the derivative ${\displaystyle \ F'(x)~.}$ Note that for arguments ${\displaystyle \ x\ }$ close to some fixed point ${\displaystyle \ x_{\star }\ ,}$ fixed point functions ${\displaystyle \ F\ }$ and their linear operators ${\displaystyle \ G\ }$ meeting condition (1), ${\displaystyle \ F'(x)\ \approx \ G{\bigl (}F\left(x\right),x{\bigr )}\ \approx \ I\ ,}$ where ${\displaystyle \ I\ }$ is the identity operator.

In the case that division is possible in the Banach space, the generalized iteration formula is given by

${\displaystyle x_{n+1}=x_{n}+{\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}^{-1}{\Bigl [}\ F\left(x_{n}\right)-x_{n}\ {\Bigr ]}\ ,}$

for ${\displaystyle \ n=1,\ 2,\ 3,\ ...~.}$ In the more general case in which division may not be possible, the iteration formula requires finding a solution ${\displaystyle \ x_{n+1}\ }$ close to ${\displaystyle \ x_{n}\ }$ for which

${\displaystyle {\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}{\bigl (}\ x_{n+1}-x_{n}\ {\bigr )}=F\left(x_{n}\right)-x_{n}~.}$

Equivalently, one may seek the solution ${\displaystyle \ x_{n+1}\ }$ to the somewhat reduced form

${\displaystyle {\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}\ x_{n+1}={\Bigl [}\ F\left(x_{n}\right)-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ x_{n}\ {\Bigr ]}\ ,}$

with all the values inside square brackets being independent of ${\displaystyle \ x_{n+1}\ :}$ The bracketed terms all only depend on ${\displaystyle \ x_{n}\ }$. However, the second form may not be as numerically stable as the first; because the first form involves finding a value for a (hopefully) small difference, it may be numerically more likely to avoid excessively large or erratic changes to the iterated value ${\displaystyle \ x_{n}~.}$

If the linear operator ${\displaystyle \ G\ }$ satisfies

${\displaystyle {\Bigl \|}G\left(u,v\right)-G\left(x,y\right){\Bigr \|}\leq k{\biggl (}{\Bigl \|}u-x{\Bigr \|}+{\Bigr \|}v-y{\Bigr \|}{\biggr )}}$

for some positive real constant ${\displaystyle \ k\ ,}$ then the method converges quadratically to a fixed point of ${\displaystyle \ F\ }$ if the initial approximation ${\displaystyle \ x_{0}\ }$ is "sufficiently close" to the desired solution ${\displaystyle \ x_{\star }\ }$ that satisfies ${\displaystyle \ x_{\star }=F(x_{\star })~.}$

Notes

1. ^ For rare special case functions the derivative for Newton's method can be calculated at negligible cost, by using saved parts from evaluation of the main function. If optimized in this way, Newton's method becomes only slightly more costly per step than the secant method, and benefits from slightly faster convergence.
2. ^ ^{a b} The condition ${\displaystyle -1<f'(x_{\star })<0\ }$ ensures that if ${\displaystyle \ f\ }$ was used as a correction-function for ${\displaystyle \ x\ ,}$ for finding its own solution, it would step in the direction of the solution (${\displaystyle \ f'<0\ }$), and that the new value would tend to lie in between the solution and the prior value (${\displaystyle -1\ <f'\ }$). But note that ${\displaystyle \ f\ }$ is only a self-correction function in principle. It is not actually used for that purpose, and it is not required to be efficient, even if it were so used.
3. ^ The divided difference ${\displaystyle \ g\ }$ is either a forward-type or backward-type
   divided difference
   , depending on the sign of ${\displaystyle \ h\ }$.
4. ^ Because ${\displaystyle \ f(x_{n}+h)\ }$ requires the prior calculation of ${\displaystyle \ h\equiv f(x_{n})\ ,}$ the two evaluations must be done sequentially – the algorithm per se cannot be made faster by running the function evaluations in parallel. This is yet another disadvantage of Steffensen's method.

References

1. ^ ^{a b c} Dahlquist, Germund; Björck, Åke (1974). Numerical Methods. Translated by Anderson, Ned. Englewood Cliffs, NJ: Prentice Hall. pp. 230–231.
2. ^ Johnson, L.W.; Scholz, D.R. (June 1968). "On Steffensen's method".
   JSTOR 2949443
   .