Steffensen's method

Steffensen's method is a numerical method to find the roots of functions. It is similar to Newton's method, but, unlike Newton's, does not require derivatives. Like Newton's, however, it is prone to diverge from finding a solution.

In this task we will do a root-finding problem that illustrates both the advantage—that no derivatives are required—and the disadvantage—that the method easily diverges. We will try to find intersection points of two Bézier curves.

Steffensen's method

We will be using the variant of Steffensen's method that employs Aitken's extrapolation. Aitken's extrapolation is illustrated by the following ATS function:

fun aitken
     (f  : double -> double,   (* function double to double *)
      p0 : double)             (* initial fixed point estimate *)
   : double =
 let
   val p1 = f(p0)
   val p2 = f(p1)
   val p1m0 = p1 - p0
 in
   p0 - (p1m0 * p1m0) / (p2 - (2.0 * p1) + p0)
 end

The return value is a function of $p0$, $p1=f(p0)$, and $p2=f(f(p0))$. What one is trying to find is a so-called fixed point of $f$: a value of $p$ such that $f(p)=p$. Our implementation of Steffensen's method will be to repeat Aitken's extrapolation until either a tolerance is satisfied or too many iterations have been executed:

fun steffensen_aitken     (* finds fixed point p such that f(p) = p *)
     (f       : double -> double, (* function double to double *)
      pinit   : double,           (* initial estimate *)
      tol     : double,           (* tolerance *)
      maxiter : int)              (* maximum number of iterations *)
   : Option (double) =     (* return a double, IF tolerance is met *)
 let
   var p0   : double = pinit
   var p    : double = aitken (f, p0)
   var iter : int = 1          (* iteration counter *)
 in
   while (abs (p - p0) > tol andalso iter < maxiter)
     begin
       p0 := p;
       p := aitken (f, p0);
       iter := iter + 1
     end;
   if abs (p - p0) > tol then None () else Some (p)
 end

The function steffensen_aitken will find a fixed point for us, but how can we use that to find a root? In particular, suppose one wants to find $t$ such that $g(t)=0$. Then what one can do (and what we will try to do) is find a fixed point $p$ of $f(t)=g(t)+t$. For then $p=f(p)=g(p)+p$, and therefore $g(p)=0$.