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$.

The problem we wish to solve

Suppose we have two quadratic planar Bézier curves, with control points $(-1,0),(0,10),(1,0)$ and $(2,1),(-8,2),(2,3)$, respectively. These two parabolas are shown in the following figure. As you can see, they intersect at four points.

We want to try to find the points of intersection.

The method we will use is implicitization. In this method, one first rewrites one of the curves as an implicit equation in $x$ and $y$. For this we will use the parabola that is convex upwards: it has implicit equation $5x^{2}+y-5=0$. Then what one does is plug the parametric equations of the other curve into the implicit equation. The resulting equation is $5(x(t))^{2}+y(t)-5=0$, where $t$ is the independent parameter for the curve that is convex leftwards. After expansion, this will be a degree-4 equation in $t$. Find its four roots and you have found the points of intersection of the two curves.

That is easier said than done.