Steffensen's method

From Rosetta Code
Task
Steffensen's method
You are encouraged to solve this task according to the task description, using any language you may know.

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 towards not finding a solution at all.

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 often gives no answer. 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 , , and . What one is trying to find is a so-called fixed point of : a value of such that . 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 such that . Then what one can do (and what we will try to do) is find a fixed point of . In that case, , and therefore .

The problem we wish to solve

Suppose we have two quadratic planar Bézier curves, with control points and , respectively. These two parabolas are shown in the following figure. As you can see, they intersect at four points.


Two parabolas in the plane that intersect at four points, plus their Bézier control polygons.


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 and . For this we will use the parabola that is convex upwards: it has implicit equation . Then what one does is plug the parametric equations of the other curve into the implicit equation. The resulting equation is , where is the independent parameter for the curve that is convex leftwards. After expansion, this will be a degree-4 equation in . Find its four roots and you have found the points of intersection of the two curves.

That is easier said than done.

There are excellent ways to solve this problem, but we will not be using those. Our purpose, after all, is to illustrate Steffensen's method: both its advantages and its drawbacks. Let us look at an advantage: to use Steffensen's method (which requires only function values, not derivatives), we do not actually have to expand that polynomial in . Instead, we can evaluate and and plug the resulting numbers into the implicit equation. What is more, we do not even need to write and as polynomials, but instead can evaluate them directly from their control points, using de Casteljau's algorithm:

fun de_casteljau
     (c0 : double,         (* control point coordinates (one axis) *)
      c1 : double,
      c2 : double,
      t  : double)             (* the independent parameter *)
   : double =                  (* value of x(t) or y(t) *)
 let
   val s = 1.0 - t
   val c01 = (s * c0) + (t * c1)
   val c12 = (s * c1) + (t * c2)
   val c012 = (s * c01) + (t * c12)
 in
   c012
 end

fun x_convex_left_parabola (t : double) : double =
 de_casteljau (2.0, ~8.0, 2.0, t)

fun y_convex_left_parabola (t : double) : double =
 de_casteljau (1.0, 2.0, 3.0, t)

Plugging and into the implicit equation, and writing as the function whose fixed points are to be found:

fun implicit_equation (x : double, y : double) : double =
 (5.0 * x * x) + y - 5.0

fun f (t : double) : double = (* find fixed points of this function *)
 let
   val x = x_convex_left_parabola (t)
   val y = y_convex_left_parabola (t)
 in
   implicit_equation (x, y) + t
 end

What is left to be done is simple, but tricky: use the steffensen_aitken function to find those fixed points. This is where a huge disadvantage of Steffensen's method (shared with Newton's method) will be illustrated. What is likely to happen is that only some of the intersection points will be found. Furthermore, for most of our initial estimates, Steffensen's method will not give us an answer at all.

The task

Use the methods described above to try to find intersection points of the two parabolas. For initial estimates, use . Choose reasonable settings for tolerance and maximum iterations. (The ATS example has 0.00000001 and 1000, respectively.) For each initial estimate, if tolerance was not met before maximum iterations was reached, print that there was no answer. On the other hand, if there was an answer, plug the resulting value of into and to find the intersection point . Check that that this answer is correct, by plugging into the implicit equation. If it is not correct, print that Steffensen's method gave a spurious answer. Otherwise print the value of .

(The ATS example has purposely been written in a fashion to be readable as if it were pseudocode. You may use it as a reference implementation.)

See also

ALGOL 68

Translation of: ATS

Uses a UNION(VOID,REAL) mode to represent the optional REAL value.

BEGIN # Steffenses's method - translated from the ATS sample #

    MODE OPTIONALREAL = UNION( VOID, REAL );

    # Aitken's extrapolation #
    PROC aitken = ( PROC(REAL)REAL f   # function double to double #
                  , REAL           p0  # initial fixed point estimate #
                  ) REAL:
         BEGIN
            REAL p1   = f( p0 );
            REAL p2   = f( p1 );
            REAL p1m0 = p1 - p0;
            p0 - ( p1m0 * p1m0 ) / ( p2 - ( 2 * p1 ) + p0 )
         END;

    # finds fixed point p such that f(p) = p #
    PROC steffensen aitken = ( PROC(REAL)REAL f # function double to double #
                             , REAL pinit       # initial estimate #
                             , REAL tol         # tolerance #
                             , INT  maxiter     # maximum number of iterations #
                             ) OPTIONALREAL:    # return a REAL, IF tolerance is met #
         BEGIN
            REAL p0  := pinit;
            REAL p   := aitken( f, p0 );        # first iteration #
            TO max iter - 1 WHILE ABS ( p - p0 ) > tol DO
                p0   := p;
                p    := aitken( f, p0 )
            OD;
            IF ABS ( p - p0 ) > tol THEN EMPTY ELSE p FI
         END;

    PROC de casteljau = ( REAL c0      # control point coordinates (one axis) #
                        , REAL c1
                        , REAL c2
                        , REAL t       # the independent parameter #
                        ) REAL:        # value of x(t) or y(t) #
         BEGIN
            REAL s   = 1 - t;
            REAL c01 = ( s * c0 ) + ( t * c1 );
            REAL c12 = ( s * c1 ) + ( t * c2 );
            ( s * c01 ) + ( t * c12 )
         END;

    PROC x convex left parabola = ( REAL t )REAL: de casteljau( 2, -8, 2, t );
    PROC y convex left parabola = ( REAL t )REAL: de casteljau( 1, 2, 3, t );

    PROC implicit equation = ( REAL x, y )REAL: ( 5 * x * x ) + y - 5;

    PROC f = ( REAL t )REAL: # find fixed points of this function #
         BEGIN
            REAL x = x convex left parabola( t );
            REAL y = y convex left parabola( t );
            implicit equation( x, y ) + t
         END;

    # returns v formatted with 1 digit before thw point and 6 after, with a leading sign only if negative #
    OP FMT = ( REAL v )STRING: fixed( v, IF v < 0 THEN -9 ELSE -8 FI, 6 );

    BEGIN
        REAL t0 := 0;
        TO 11 DO
            print( ( "t0 = ", FMT t0, " : " ) );
            CASE steffensen aitken( f, t0, 0.00000001, 1000 )
              IN ( VOID   ): print( ( "no answer", newline) )
               , ( REAL t ):
                     IF REAL x = x convex left parabola( t )
                           , y = y convex left parabola( t )
                           ;
                        ABS implicit equation( x, y ) <= 0.000001
                     THEN
                         print( ( "intersection at (", FMT x, ", ", FMT y, ")", newline ) )
                     ELSE
                         print( ( "spurious solution", newline ) )
                     FI
            ESAC;
            t0 +:= 0.1
        OD
    END

END
Output:

Same as the ATS sample:

t0 = 0.000000 : no answer
t0 = 0.100000 : intersection at (0.881025, 1.118975)
t0 = 0.200000 : no answer
t0 = 0.300000 : no answer
t0 = 0.400000 : no answer
t0 = 0.500000 : no answer
t0 = 0.600000 : no answer
t0 = 0.700000 : no answer
t0 = 0.800000 : no answer
t0 = 0.900000 : intersection at (-0.681025, 2.681025)
t0 = 1.000000 : no answer

ATS

#include "share/atspre_staload.hats"

fun aitken                      (* Aitken's extrapolation *)
      (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

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

fun de_casteljau
      (c0 : double,         (* control point coordinates (one axis) *)
       c1 : double,
       c2 : double,
       t  : double)             (* the independent parameter *)
    : double =                  (* value of x(t) or y(t) *)
  let
    val s = 1.0 - t
    val c01 = (s * c0) + (t * c1)
    val c12 = (s * c1) + (t * c2)
    val c012 = (s * c01) + (t * c12)
  in
    c012
  end

fun x_convex_left_parabola (t : double) : double =
  de_casteljau (2.0, ~8.0, 2.0, t)

fun y_convex_left_parabola (t : double) : double =
  de_casteljau (1.0, 2.0, 3.0, t)

fun implicit_equation (x : double, y : double) : double =
  (5.0 * x * x) + y - 5.0

fun f (t : double) : double = (* find fixed points of this function *)
  let
    val x = x_convex_left_parabola (t)
    val y = y_convex_left_parabola (t)
  in
    implicit_equation (x, y) + t
  end

implement main0 () =
  let
    var i : int = 0
    var t0 : double = 0.0
  in
    while (not (i = 11))
      begin
        begin
          print! ("t0 = ", t0, " : ");
          case steffensen_aitken (f, t0, 0.00000001, 1000) of
          | None () => println! ("no answer")
          | Some (t) =>
            let
              val x = x_convex_left_parabola (t)
              val y = y_convex_left_parabola (t)
            in
              if abs (implicit_equation (x, y)) <= 0.000001 then
                println! ("intersection at (", x, ", ", y, ")")
              else
                (* In my experience, it is possible for the algorithm
                   to achieve tolerance and yet give a spurious
                   answer. Exploration of this phenomenon is beyond
                   the scope of this task. Such spurious answers are
                   easy to discard. *)
                println! ("spurious solution")
            end
        end;
        i := i + 1;
        t0 := t0 + 0.1
      end
  end
Output:
t0 = 0.000000 : no answer
t0 = 0.100000 : intersection at (0.881025, 1.118975)
t0 = 0.200000 : no answer
t0 = 0.300000 : no answer
t0 = 0.400000 : no answer
t0 = 0.500000 : no answer
t0 = 0.600000 : no answer
t0 = 0.700000 : no answer
t0 = 0.800000 : no answer
t0 = 0.900000 : intersection at (-0.681025, 2.681025)
t0 = 1.000000 : no answer

C

Translation of: ATS
#include <stdio.h>
#include <math.h>

double aitken(double (*f)(double), double p0) {
    double p1 = f(p0);
    double p2 = f(p1);
    double p1m0 = p1 - p0;
    return p0 - p1m0 * p1m0 / (p2 - 2.0 * p1 + p0);
}

double steffensenAitken(double (*f)(double), double pinit, double tol, int maxiter) {
    double p0 = pinit;
    double p = aitken(f, p0);
    int iter = 1;
    while (fabs(p - p0) > tol && iter < maxiter) {
        p0 = p;
        p = aitken(f, p0);
        ++iter;
    }
    if (fabs(p - p0) > tol) return nan("");
    return p;
}

double deCasteljau(double c0, double c1, double c2, double t) {
    double s = 1.0 - t;
    double c01 = s * c0 + t * c1;
    double c12 = s * c1 + t * c2;
    return s * c01 + t * c12;
}

double xConvexLeftParabola(double t) {
    return deCasteljau(2.0, -8.0, 2.0, t);
}

double yConvexRightParabola(double t) {
    return deCasteljau(1.0, 2.0, 3.0, t);
}

double implicitEquation(double x, double y) {
    return 5.0 * x * x + y - 5.0;
}

double f(double t) {
    double x = xConvexLeftParabola(t);
    double y = yConvexRightParabola(t);
    return implicitEquation(x, y) + t;
}

int main() {
    double t0 = 0.0, t, x, y;
    int i;
    for (i = 0; i < 11; ++i) {
        printf("t0 = %0.1f : ", t0);
        t = steffensenAitken(f, t0, 0.00000001, 1000);
        if (isnan(t)) {
            printf("no answer\n");
        } else {
            x = xConvexLeftParabola(t);
            y = yConvexRightParabola(t);
            if (fabs(implicitEquation(x, y)) <= 0.000001) {
                printf("intersection at (%f, %f)\n", x, y);
            } else {
                printf("spurious solution\n");
            }
        }
        t0 += 0.1;
    }
    return 0;
}
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at (0.881025, 1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.681025, 2.681025)
t0 = 1.0 : no answer

C++

Translation of: C
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

double aitken(double (*f)(double), double p0) {
    double p1 = f(p0);
    double p2 = f(p1);
    double p1m0 = p1 - p0;
    return p0 - p1m0 * p1m0 / (p2 - 2.0 * p1 + p0);
}

double steffensenAitken(double (*f)(double), double pinit, double tol, int maxiter) {
    double p0 = pinit;
    double p = aitken(f, p0);
    int iter = 1;
    while (fabs(p - p0) > tol && iter < maxiter) {
        p0 = p;
        p = aitken(f, p0);
        ++iter;
    }
    if (fabs(p - p0) > tol) return nan("");
    return p;
}

double deCasteljau(double c0, double c1, double c2, double t) {
    double s = 1.0 - t;
    double c01 = s * c0 + t * c1;
    double c12 = s * c1 + t * c2;
    return s * c01 + t * c12;
}

double xConvexLeftParabola(double t) {
    return deCasteljau(2.0, -8.0, 2.0, t);
}

double yConvexRightParabola(double t) {
    return deCasteljau(1.0, 2.0, 3.0, t);
}

double implicitEquation(double x, double y) {
    return 5.0 * x * x + y - 5.0;
}

double f(double t) {
    double x = xConvexLeftParabola(t);
    double y = yConvexRightParabola(t);
    return implicitEquation(x, y) + t;
}

int main() {
    double t0 = 0.0, t, x, y;
    int i;
    for (i = 0; i < 11; ++i) {
        cout << "t0 = " << setprecision(1) << t0 << " : ";
        t = steffensenAitken(f, t0, 0.00000001, 1000);
        if (isnan(t)) {
            cout << "no answer << endl;
        } else {
            x = xConvexLeftParabola(t);
            y = yConvexRightParabola(t);
            if (fabs(implicitEquation(x, y)) <= 0.000001) {
                cout << "intersection at (" << x << ", " << y << ")" << endl;
            } else {
                cout << "spurious solution" << endl;
            }
        }
        t0 += 0.1;
    }
    return 0;
}
Output:
Same as C example.

EasyLang

Translation of: C
func deCasteljau c0 c1 c2 t .
   s = 1 - t
   c01 = s * c0 + t * c1
   c12 = s * c1 + t * c2
   return s * c01 + t * c12
.
func xConvexLeftPar t .
   return deCasteljau 2 (-8) 2 t
.
func yConvexRightPar t .
   return deCasteljau 1 2 3 t
.
func implicitEq x y .
   return 5 * x * x + y - 5
.
func f t r .
   x = xConvexLeftPar t
   y = yConvexRightPar t
   r = implicitEq x y
   return r + t
.
func aitken p0 .
   p1 = f p0 p1
   p2 = f p1 p2
   p1m0 = p1 - p0
   return p0 - p1m0 * p1m0 / (p2 - 2 * p1 + p0)
.
func steffAitken p0 tol maxiter .
   for i to maxiter
      p = aitken p0
      if abs (p - p0) < tol
         return p
      .
      p0 = p
   .
   return number "nan"
.
for i to 11
   numfmt 1 0
   write "t0 = " & t0 & " : "
   t = steffAitken t0 0.00000001 1000
   numfmt 3 0
   if t <> t
      # nan
      print "no answer"
   else
      x = xConvexLeftPar t
      y = yConvexRightPar t
      r = implicitEq x y
      if abs r <= 0.000001
         print "intersection at (" & x & " " & y & ")"
      else
         print "spurious solution"
      .
   .
   t0 += 0.1
.

FreeBASIC

Translation of: ATS
Function aitken(f As Function(As Double) As Double, p0 As Double) As Double
    Dim As Double p1 = f(p0)
    Dim As Double p2 = f(p1)
    Dim As Double p1m0 = p1 - p0    
    
    Return p0 - (p1m0 * p1m0) / (p2 - (2 * p1) + p0)
End Function

Function steffensenAitken(f As Function(As Double) As Double, pinit As Double, tol As Double, maxiter As Integer) As Double
    Dim As Double p0 = pinit
    Dim As Double p = aitken(f, p0)
    Dim As Integer iter = 1
    While Abs(p-p0) > tol Andalso iter < maxiter
        p0 = p
        p = aitken (f, p0)
        iter += 1
    Wend
    
    If Abs (p - p0) > tol Then Return 0 Else Return p
End Function

Function deCasteljau(c0 As Double, c1 As Double, c2 As Double, t As Double) As Double
    Dim As Double s = 1 - t
    Dim As Double c01 = (s * c0) + (t * c1)
    Dim As Double c12 = (s * c1) + (t * c2)
    Dim As Double c012 = (s * c01) + (t * c12)    
    Return c012
End Function

Function xConvexLeftParabola(t As Double) As Double
    Return deCasteljau(2, -8, 2, t)
End Function

Function yConvexLeftParabola(t As Double) As Double
    Return deCasteljau(1, 2, 3, t)
End Function

Function implicitEquation(x As Double, y As Double) As Double 
    Return (5 * x * x) + y - 5
End Function

Function f(t As Double) As Double
    Dim As Double x = xConvexLeftParabola(t)
    Dim As Double y = yConvexLeftParabola(t)
    
    Return implicitEquation(x, y) + t
End Function

Dim As Double t0 = 0

For i As Integer = 0 To 10
    Print Using "t0 = #.# : "; t0;
    Dim As Double t = steffensenAitken(@f, t0, 0.00000001, 1000)
    If t = 0 Then
        Print "no answer"
    Else
        Dim As Double x = xConvexLeftParabola(t)
        Dim As Double y = yConvexLeftParabola(t)
        
        If Abs(implicitEquation(x, y)) <= 0.000001 Then
            Print Using "intersection at (##.######, ##.######)"; x; y
        Else
            Print "spurious solution"
        End If
    End If
    t0 += 0.1
Next

Sleep
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at ( 0.881025,  1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.681025,  2.681025)
t0 = 1.0 : no answer

Go

Translation of: C
package main

import (
    "fmt"
    "math"
)

type fun = func(float64) float64

func aitken(f fun, p0 float64) float64 {
    p1 := f(p0)
    p2 := f(p1)
    p1m0 := p1 - p0
    return p0 - p1m0*p1m0/(p2-2.0*p1+p0)
}

func steffensenAitken(f fun, pinit, tol float64, maxiter int) float64 {
    p0 := pinit
    p := aitken(f, p0)
    iter := 1
    for math.Abs(p-p0) > tol && iter < maxiter {
        p0 = p
        p = aitken(f, p0)
        iter++
    }
    if math.Abs(p-p0) > tol {
        return math.NaN()
    }
    return p
}

func deCasteljau(c0, c1, c2, t float64) float64 {
    s := 1.0 - t
    c01 := s*c0 + t*c1
    c12 := s*c1 + t*c2
    return s*c01 + t*c12
}

func xConvexLeftParabola(t float64) float64 {
    return deCasteljau(2.0, -8.0, 2.0, t)
}

func yConvexRightParabola(t float64) float64 {
    return deCasteljau(1.0, 2.0, 3.0, t)
}

func implicitEquation(x, y float64) float64 {
    return 5.0*x*x + y - 5.0
}

func f(t float64) float64 {
    x := xConvexLeftParabola(t)
    y := yConvexRightParabola(t)
    return implicitEquation(x, y) + t
}

func main() {
    t0 := 0.0
    for i := 0; i < 11; i++ {
        fmt.Printf("t0 = %0.1f : ", t0)
        t := steffensenAitken(f, t0, 0.00000001, 1000)
        if math.IsNaN(t) {
            fmt.Println("no answer")
        } else {
            x := xConvexLeftParabola(t)
            y := yConvexRightParabola(t)
            if math.Abs(implicitEquation(x, y)) <= 0.000001 {
                fmt.Printf("intersection at (%f, %f)\n", x, y)
            } else {
                fmt.Println("spurious solution")
            }
        }
        t0 += 0.1
    }
}
Output:
Same as C example

Java

Translation of: ATS
import java.util.Optional;

public class Steffensen {
    static double aitken(double p0) {
        double p1 = f(p0);
        double p2 = f(p1);
        double p1m0 = p1 - p0;
        return p0 - p1m0 * p1m0 / (p2 - 2.0 * p1 + p0);
    }

    static Optional<Double> steffensenAitken(double pinit, double tol, int maxiter) {
        double p0 = pinit;
        double p = aitken(p0);
        int iter = 1;
        while (Math.abs(p - p0) > tol && iter < maxiter) {
            p0 = p;
            p = aitken(p0);
            iter++;
        }
        if (Math.abs(p - p0) > tol) return Optional.empty();
        return Optional.of(p);
    }

    static double deCasteljau(double c0, double c1, double c2, double t) {
        double s = 1.0 - t;
        double c01 = s * c0 + t * c1;
        double c12 = s * c1 + t * c2;
        return s * c01 + t * c12;
    }

    static double xConvexLeftParabola(double t) {
        return deCasteljau(2.0, -8.0, 2.0, t);
    }

    static double yConvexRightParabola(double t) {
        return deCasteljau(1.0, 2.0, 3.0, t);
    }

    static double implicitEquation(double x, double y) {
        return 5.0 * x * x + y - 5.0;
    }

    static double f(double t) {
        double x = xConvexLeftParabola(t);
        double y = yConvexRightParabola(t);
        return implicitEquation(x, y) + t;
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        for (int i = 0; i < 11; ++i) {
            System.out.printf("t0 = %3.1f : ", t0);
            Optional<Double> t = steffensenAitken(t0, 0.00000001, 1000);
            if (!t.isPresent()) {
                System.out.println("no answer");
            } else {
                double x = xConvexLeftParabola(t.get());
                double y = yConvexRightParabola(t.get());
                if (Math.abs(implicitEquation(x, y)) <= 0.000001) {
                    System.out.printf("intersection at (%f, %f)\n", x, y);
                } else {
                    System.out.println("spurious solution");
                }
            }
            t0 += 0.1;
        }
    }
}
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at (0.881025, 1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.681025, 2.681025)
t0 = 1.0 : no answer

jq

Adapted from Wren

Works with: jq

Also works with gojq, the Go implementation of jq

def aitken(f):
  . as $p0
  | f as $p1
  | ($p1|f) as $p2
  | ($p1 - $p0) as $p1m0
  | $p0 - $p1m0 * $p1m0 / ($p2 - 2 * $p1 + $p0);

def steffensenAitken(f; $pinit; $tol; $maxiter):
  { p0:  $pinit}
  | .p = (.p0|aitken(f))
  | .iter = 1
  | until( ((.p - .p0)|length) <= $tol or .iter >= $maxiter;
      .p0 = .p
      | .p = (.p0|aitken(f))
      | .iter += 1 )
  | if ((.p - .p0)|length > $tol) then null else .p end;

def deCasteljau($c0; $c1; $c2):
  (1 - .) as $s
  | ($s * $c0 + . * $c1) as $c01
  | ($s * $c1 + . * $c2) as $c12
  | $s * $c01 + . * $c12;

def xConvexLeftParabola:  deCasteljau(2; -8; 2);
def yConvexRightParabola: deCasteljau(1;  2; 3);

def implicitEquation($x; $y): 5 * $x * $x + $y - 5;

def f:
    xConvexLeftParabola as $x
  | yConvexRightParabola as $y
  | implicitEquation($x; $y) + . ;

def pp: . * 100 | round / 100;

foreach range (0; 11) as $i ( {t0:0};
  .emit = "t0 = \(.t0|pp): "
  | steffensenAitken(f; .t0; 0.00000001; 1000) as $t
  | if $t
    then ($t | xConvexLeftParabola) as $x
    |    ($t | yConvexRightParabola) as $y
    | if (implicitEquation($x; $y)|length) <= 0.000001
      then .emit += "intersection at \($x), \($y)"
      else .emit += "spurious solution"
      end
    else .emit += "no answer"
    end
  | .t0 += 0.1 )
| .emit
Output:
t0 = 0: no answer
t0 = 0.1: intersection at 0.8810249675906652, 1.1189750324093346
t0 = 0.2: no answer
t0 = 0.3: no answer
t0 = 0.4: no answer
t0 = 0.5: no answer
t0 = 0.6: no answer
t0 = 0.7: no answer
t0 = 0.8: no answer
t0 = 0.9: intersection at -0.6810249675905098, 2.6810249675906883
t0 = 1: no answer

Julia

Translation of: C
""" Aitken's extrapolation """
function aitken(f, p0)
    p1 = f(p0)
    p2 = f(p1)
    return p0 - (p1 - p0)^2 / (p2 - 2 * p1 + p0)
end

""" Steffensen's method using Aitken """
function steffensen_aitken(f, pinit, tol, maxiter)
    p0 = pinit
    p = aitken(f, p0)
    iter = 1
    while abs(p - p0) > tol && iter < maxiter
        p0 = p
        p = aitken(f, p0)
        iter += 1
    end
    return abs(p - p0) > tol ? NaN : p
end

""" deCasteljau function """
function deCasteljau(c0, c1, c2, t)
    s = 1.0 - t
    return s * (s * c0 + t * c1) + t * (s * c1 + t * c2)
end

xConvexLeftParabola(t) = deCasteljau(2, -8, 2, t)
yConvexRightParabola(t) = deCasteljau(1, 2, 3, t)
implicit_equation(x, y) = 5 * x^2 + y - 5

""" may return NaN on overflow """
function f(t)
    return t in [nothing, NaN, Inf, -Inf] ? NaN :
       implicit_equation(xConvexLeftParabola(t), yConvexRightParabola(t)) + t
end

""" test the example """
function test_steffensen(tol = 0.00000001, iters = 1000, stepsize = 0.1)
    for t0 in 0:stepsize:1.1
        print("t0 = $t0 : ")
        t = steffensen_aitken(f, t0, tol, iters)
        if isnan(t)
            println("no answer")
        else
            x = xConvexLeftParabola(t)
            y = yConvexRightParabola(t)
            if abs(implicit_equation(x, y)) <= tol
                println("intersection at ($(Float32(x)), $(Float32(y)))")
            else
                println("spurious solution")
            end
        end
    end
    return 0
end

test_steffensen()
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at (0.88102496, 1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.68102497, 2.681025)
t0 = 1.0 : no answer
t0 = 1.1 : no answer

MATLAB

Translation of: Julia
clear all;close all;clc;
tol = 0.00000001; iters = 1000; stepsize = 0.1;
test_steffensen(tol, iters, stepsize);


% Aitken's extrapolation
function p = aitken(f, p0)
    p1 = f(p0);
    p2 = f(p1);
    p = p0 - (p1 - p0)^2 / (p2 - 2 * p1 + p0);
end

% Steffensen's method using Aitken
function p = steffensen_aitken(f, pinit, tol, maxiter)
    p0 = pinit;
    p = aitken(f, p0);
    iter = 1;
    while abs(p - p0) > tol && iter < maxiter
        p0 = p;
        p = aitken(f, p0);
        iter = iter + 1;
    end
    if abs(p - p0) > tol
        p = NaN;
    end
end

% deCasteljau function
function result = deCasteljau(c0, c1, c2, t)
    s = 1.0 - t;
    result = s * (s * c0 + t * c1) + t * (s * c1 + t * c2);
end

% Helper functions
function result = xConvexLeftParabola(t)
    result = deCasteljau(2, -8, 2, t);
end

function result = yConvexRightParabola(t)
    result = deCasteljau(1, 2, 3, t);
end

function result = implicit_equation(x, y)
    result = 5 * x^2 + y - 5;
end

% Main function
function result = f(t)
    if isnan(t) || isinf(t)
        result = NaN;
    else
        result = implicit_equation(xConvexLeftParabola(t), yConvexRightParabola(t)) + t;
    end
end

% Test example
function test_steffensen(tol, iters, stepsize)
    if nargin < 1
        tol = 0.00000001;
    end
    if nargin < 2
        iters = 1000;
    end
    if nargin < 3
        stepsize = 0.1;
    end

    for t0 = 0:stepsize:1.1
        fprintf('t0 = %f : ', t0);
        t = steffensen_aitken(@f, t0, tol, iters);
        if isnan(t)
            fprintf('no answer\n');
        else
            x = xConvexLeftParabola(t);
            y = yConvexRightParabola(t);
            if abs(implicit_equation(x, y)) <= tol
                fprintf('intersection at (%f, %f)\n', x, y);
            else
                fprintf('spurious solution\n');
            end
        end
    end
end
Output:
t0 = 0.000000 : no answer
t0 = 0.100000 : intersection at (0.881025, 1.118975)
t0 = 0.200000 : no answer
t0 = 0.300000 : no answer
t0 = 0.400000 : no answer
t0 = 0.500000 : no answer
t0 = 0.600000 : no answer
t0 = 0.700000 : no answer
t0 = 0.800000 : no answer
t0 = 0.900000 : intersection at (-0.681025, 2.681025)
t0 = 1.000000 : no answer
t0 = 1.100000 : no answer

Nim

Translation of: C
import std/[math, strformat]

proc aitken(f: proc(x: float): float; p0: float): float =
  let p1 = f(p0)
  let p2 = f(p1)
  let p1m0 = p1 - p0
  result = p0 - p1m0 * p1m0 / (p2 - 2 * p1 + p0)

proc steffensenAitken(f: proc(x: float): float; pinit, tol: float; maxiter: int): float =
  var p0 = pinit
  result = aitken(f, p0)
  var iter = 1
  while abs(result - p0) > tol and iter < maxiter:
    p0 = result
    result = aitken(f, p0)
    inc iter
  if abs(result - p0) > tol: return Nan

func deCasteljau(c0, c1, c2, t: float): float =
  let s = 1 - t
  let c01 = s * c0 + t * c1
  let c12 = s * c1 + t * c2
  result = s * c01 + t * c12

template xConvexLeftParabola(t: float): float =
  deCasteljau(2, -8, 2, t)

template yConvexRightParabola(t: float): float =
  deCasteljau(1, 2, 3, t)

func implicitEquation(x, y: float): float =
  5 * x * x + y - 5

func f(t: float): float =
  let x = xConvexLeftParabola(t)
  let y = yConvexRightParabola(t)
  result = implicitEquation(x, y) + t

var t0 = 0.0
var x, y: float
for i in 0..10:
  stdout.write &"t0 = {t0:0.1f} → "
  let t = steffensenAitken(f, t0, 0.00000001, 1000)
  if t.isNan:
    echo "no answer"
  else:
    x = xConvexLeftParabola(t);
    y = yConvexRightParabola(t);
    if abs(implicitEquation(x, y)) <= 0.000001:
      echo &"intersection at ({x:.6f}, {y:.6f})"
    else:
      echo "spurious solution"
  t0 += 0.1
Output:
t0 = 0.0 → no answer
t0 = 0.1 → intersection at (0.881025, 1.118975)
t0 = 0.2 → no answer
t0 = 0.3 → no answer
t0 = 0.4 → no answer
t0 = 0.5 → no answer
t0 = 0.6 → no answer
t0 = 0.7 → no answer
t0 = 0.8 → no answer
t0 = 0.9 → intersection at (-0.681025, 2.681025)
t0 = 1.0 → no answer

Phix

Translation of: C
with javascript_semantics
function aitken(integer f, atom p0)
    atom p1 = f(p0),
         p2 = f(p1),
       p1m0 = p1 - p0   
    return p0 - (p1m0 * p1m0) / (p2 - (2 * p1) + p0)
end function

function steffensenAitken(integer f, maxiter, atom pinit, tol)
    atom p0 = pinit, p = aitken(f, p0)
    integer iter = 1
    while abs(p-p0)>tol and iter<maxiter do
        p0 = p
        p = aitken(f, p0)
        iter += 1
    end while
    return iff(abs(p-p0)>tol?"none":p)
end function

function deCasteljau(atom c0, c1, c2, t)
    atom s = 1 - t,
       c01 = (s * c0) + (t * c1),
       c12 = (s * c1) + (t * c2),
      c012 = (s * c01) + (t * c12)    
    return c012
end function

function xConvexLeftParabola(atom t)
    return deCasteljau(2, -8, 2, t)
end function

function yConvexLeftParabola(atom t)
    return deCasteljau(1, 2, 3, t)
end function

function implicitEquation(atom x, y)
    return (5 * x * x) + y - 5
end function

function f(atom t)
    atom x = xConvexLeftParabola(t),
         y = yConvexLeftParabola(t)
    return implicitEquation(x, y) + t
end function

for i=0 to 10 do
    printf(1,"t0 = %.1f : ",i/10)
    object t = steffensenAitken(f, 1000, i/10, 0.00000001)
    if string(t) then
        printf(1,"no answer\n")
    else
        atom x = xConvexLeftParabola(t),
             y = yConvexLeftParabola(t)
        if abs(implicitEquation(x, y)) <= 0.000001 then
            printf(1,"intersection at (%.6f, %.6f)\n",{x,y})
        else
            printf(1,"spurious solution\n")
        end if
    end if
end for

Output same as C

Python

Translation of: C
from math import nan, isnan
from numpy import arange


def aitken(f, p0):
    """ Aitken's extrapolation """
    p1 = f(p0)
    p2 = f(p1)
    return p0 - (p1 - p0)**2 / (p2 - 2 * p1 + p0)

def steffensen_aitken(f, pinit, tol, maxiter):
    """ Steffensen's method using Aitken """
    p0 = pinit
    p = aitken(f, p0)
    iter = 1
    while abs(p - p0) > tol and iter < maxiter:
        p0 = p
        p = aitken(f, p0)
        iter += 1
    return nan if abs(p - p0) > tol else p

def deCasteljau(c0, c1, c2, t):
    """ deCasteljau function """
    s = 1.0 - t
    return s * (s * c0 + t * c1) + t * (s * c1 + t * c2)

def xConvexLeftParabola(t): return deCasteljau(2, -8, 2, t)
def yConvexRightParabola(t): return deCasteljau(1, 2, 3, t)
def implicit_equation(x, y): return 5 * x**2 + y - 5

def f(t):
    """ Outside of NumPy arithmetic may return NoneType on overflow """
    if type(t) == type(None):
        return nan
    return implicit_equation(xConvexLeftParabola(t), yConvexRightParabola(t)) + t

def test_steffensen(tol=0.00000001, iters=1000, stepsize=0.1):
    """ test the example """
    for t0 in arange(0, 1.1, stepsize):
        print(f't0 = {t0:0.1f} : ', end='')
        t = steffensen_aitken(f, t0, tol, iters)
        if isnan(t):
            print('no answer')
        else:
            x = xConvexLeftParabola(t)
            y = yConvexRightParabola(t)
            if abs(implicit_equation(x, y)) <= tol:
                print(f'intersection at ({x}, {y})')
            else:
                print('spurious solution')
    return 0


if __name__ == '__main__':

    test_steffensen()
Output:
Output is the same as C.

Raku

Translation of: Go
# 20230928 Raku programming solution

sub aitken($f, $p0) {
   my $p2   = $f( my $p1 = $f($p0) ); 
   my $p1m0 = $p1 - $p0;
   return $p0 - $p1m0*$p1m0/($p2-2.0*$p1+$p0);
}

sub steffensenAitken($f, $pinit, $tol, $maxiter) {
   my ($iter, $p) = 1, aitken($f, my $p0 = $pinit);
   while abs($p-$p0) > $tol and $iter < $maxiter {
      $p = aitken($f, $p0 = $p);
      $iter++
   }
   return abs($p-$p0) > $tol ?? NaN !! $p
}

sub deCasteljau($c0, $c1, $c2, $t) {
   my $s = 1.0 - $t;
   return $s*($s*$c0 + $t*$c1) + $t*($s*$c1 + $t*$c2)
}

sub xConvexLeftParabola($t) { return deCasteljau(2.0, -8.0, 2.0, $t) }

sub yConvexRightParabola($t) { return deCasteljau(1.0, 2.0, 3.0, $t) }

sub implicitEquation($x, $y) { return 5.0*$x*$x + $y - 5.0 }

sub f($t) {
   implicitEquation(xConvexLeftParabola($t), yConvexRightParabola($t)) + $t
}

my $t0 = 0.0;
for ^11 {
   print "t0 = {$t0.fmt: '%0.1f'} : ";
   my $t = steffensenAitken(&f, $t0, 0.00000001, 1000);
   if $t.isNaN {
      say "no answer";
   } else {
      my ($x, $y) = xConvexLeftParabola($t), yConvexRightParabola($t);
      if abs(implicitEquation($x, $y)) <= 0.000001 {
         printf "intersection at (%f, %f)\n", $x, $y;
      } else {
         say "spurious solution";
      }
   }
   $t0 += 0.1;
}

You may Attempt This Online!

RATFOR

Translation of: ATS
function aitken (f, p0)
  implicit none

  double precision f, p0, aitken
  external f

  double precision p1, p2

  p1 = f(p0)
  p2 = f(p1)
  aitken = p0 - (p1 - p0)**2 / (p2 - (2.0d0 * p1) + p0)
end

subroutine steff (f, pinit, tol, mxiter, ierr, pfinal)
  implicit none

  double precision f, pinit, tol, pfinal
  integer mxiter, ierr
  external f

  double precision p0, p, aitken
  integer iter
  external aitken

  p0 = pinit
  p = aitken (f, p0)
  iter = 1
  while (abs (p - p0) > tol && iter < mxiter)
    {
      p0 = p
      p = aitken (f, p0)
      iter = iter + 1
    }
  if (abs (p - p0) > tol)
    ierr = -1
  else
    {
      ierr = 0
      pfinal = p
    }
end

function dcstlj (c0, c1, c2, t)
  implicit none

  double precision c0, c1, c2, t, dcstlj

  double precision s, c01, c12, c012

  s = 1.0d0 - t
  c01 = (s * c0) + (t * c1)
  c12 = (s * c1) + (t * c2)
  c012 = (s * c01) + (t * c12)
  dcstlj = c012
end

function x (t)
  implicit none
  double precision x, t, dcstlj
  external dcstlj
  x = dcstlj (2.0d0, -8.0d0, 2.0d0, t)
end

function y (t)
  implicit none
  double precision y, t, dcstlj
  external dcstlj
  y = dcstlj (1.0d0, 2.0d0, 3.0d0, t)
end

function impleq (x, y)
  implicit none
  double precision x, y, impleq
  impleq = (5.0d0 * x * x) + y - 5.0d0
end

function f (t)
  implicit none
  double precision f, t, x, y, impleq
  external x, y, impleq
  f = impleq (x(t), y(t)) + t
end

program RCstef
  implicit none

  double precision x, y, impleq, f
  external x, y, impleq, f, steff

  double precision t0, t
  integer i, ierr

  t0 = 0.0d0
  for (i = 0; i != 11; i = i + 1)
    {
      call steff (f, t0, 0.00000001d0, 1000, ierr, t)
      if (ierr < 0)
        write (*,*) "t0 = ", t0, " : no answer"
      else if (abs (impleq (x(t), y(t))) <= 0.000001)
        write (*,*) "t0 = ", t0, " : intersection at (", _
                    x(t), ", ", y(t), ")"
      else
        write (*,*) "t0 = ", t0, " : spurious solution"
      t0 = t0 + 0.1d0
    }
end
Output:

Compiled with f2c (which accounts for the ink-saving numerals):

 t0 =   0. : no answer
 t0 =   .1 : intersection at (  .881024968,   1.11897503)
 t0 =   .2 : no answer
 t0 =   .3 : no answer
 t0 =   .4 : no answer
 t0 =   .5 : no answer
 t0 =   .6 : no answer
 t0 =   .7 : no answer
 t0 =   .8 : no answer
 t0 =   .9 : intersection at ( -.681024968,   2.68102497)
 t0 =   1. : no answer

Scala

Translation of: Java
object Steffensen {

  def aitken(p0: Double): Double = {
    val p1 = f(p0)
    val p2 = f(p1)
    val p1m0 = p1 - p0
    p0 - p1m0 * p1m0 / (p2 - 2.0 * p1 + p0)
  }

  def steffensenAitken(pinit: Double, tol: Double, maxiter: Int): Option[Double] = {
    var p0 = pinit
    var p = aitken(p0)
    var iter = 1
    while (math.abs(p - p0) > tol && iter < maxiter) {
      p0 = p
      p = aitken(p0)
      iter += 1
    }
    if (math.abs(p - p0) > tol) None else Some(p)
  }

  def deCasteljau(c0: Double, c1: Double, c2: Double, t: Double): Double = {
    val s = 1.0 - t
    val c01 = s * c0 + t * c1
    val c12 = s * c1 + t * c2
    s * c01 + t * c12
  }

  def xConvexLeftParabola(t: Double): Double = deCasteljau(2.0, -8.0, 2.0, t)

  def yConvexRightParabola(t: Double): Double = deCasteljau(1.0, 2.0, 3.0, t)

  def implicitEquation(x: Double, y: Double): Double = 5.0 * x * x + y - 5.0

  def f(t: Double): Double = {
    val x = xConvexLeftParabola(t)
    val y = yConvexRightParabola(t)
    implicitEquation(x, y) + t
  }

  def main(args: Array[String]): Unit = {
    var t0 = 0.0
    for (i <- 0 until 11) {
      print(f"t0 = $t0%3.1f : ")
      steffensenAitken(t0, 0.00000001, 1000) match {
        case None => println("no answer")
        case Some(t) =>
          val x = xConvexLeftParabola(t)
          val y = yConvexRightParabola(t)
          if (math.abs(implicitEquation(x, y)) <= 0.000001) {
            println(f"intersection at ($x%f, $y%f)")
          } else {
            println("spurious solution")
          }
      }
      t0 += 0.1
    }
  }
}
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at (0.881025, 1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.681025, 2.681025)
t0 = 1.0 : no answer

Wren

Translation of: ATS
Library: Wren-fmt

Wren's only built-in numeric type is a 'double' (in C) so no surprise that we're getting the same answers as ATS when using the same tolerances and maximum number of iterations.

import "./fmt" for Fmt

var aitken = Fn.new { |f, p0|
    var p1 = f.call(p0)
    var p2 = f.call(p1)
    var p1m0 = p1 - p0
    return p0 - p1m0 * p1m0 / (p2 - 2 * p1 + p0)
}

var steffensenAitken = Fn.new { |f, pinit, tol, maxiter|
    var p0 = pinit
    var p = aitken.call(f, p0)
    var iter = 1
    while ((p - p0).abs > tol && iter < maxiter) {
        p0 = p
        p = aitken.call(f, p0)
        iter = iter + 1
    }
    if ((p - p0).abs > tol) return null
    return p
}

var deCasteljau = Fn.new { |c0, c1, c2, t|
    var s = 1 - t
    var c01 = s * c0 + t * c1
    var c12 = s * c1 + t * c2
    return s * c01 + t * c12
}

var xConvexLeftParabola  = Fn.new { |t| deCasteljau.call(2, -8, 2, t) }
var yConvexRightParabola = Fn.new { |t| deCasteljau.call(1,  2, 3, t) }

var implicitEquation = Fn.new { |x, y|  5 * x * x + y - 5 }

var f = Fn.new { |t|
    var x = xConvexLeftParabola.call(t)
    var y = yConvexRightParabola.call(t)
    return implicitEquation.call(x, y) + t
}

var t0 = 0
for (i in 0..10) {
    Fmt.write("t0 = $0.1f : ", t0)
    var t = steffensenAitken.call(f, t0, 0.00000001, 1000)
    if (!t) {
        Fmt.print("no answer")
    } else {
        var x = xConvexLeftParabola.call(t)
        var y = yConvexRightParabola.call(t)
        if (implicitEquation.call(x, y).abs <= 0.000001) {
            Fmt.print("intersection at ($f, $f)", x, y)
        } else {
            Fmt.print("spurious solution")
        }
    }
    t0 = t0 + 0.1
}
Output:
t0 = 0.0 : no answer
t0 = 0.1 : intersection at (0.881025, 1.118975)
t0 = 0.2 : no answer
t0 = 0.3 : no answer
t0 = 0.4 : no answer
t0 = 0.5 : no answer
t0 = 0.6 : no answer
t0 = 0.7 : no answer
t0 = 0.8 : no answer
t0 = 0.9 : intersection at (-0.681025, 2.681025)
t0 = 1.0 : no answer

XPL0

Translation of: C
include xpllib; \for Print
def NaN = 1e308;

func real DeCasteljau(C0, C1, C2, T);
real C0, C1, C2, T;
real S, C01, C12;
[S:= 1.0 - T;
C01:= S*C0 + T*C1;
C12:= S*C1 + T*C2;
return S*C01 + T*C12;
];

func real XConvexLeftParabola(T);
real T;
return DeCasteljau(2.0, -8.0, 2.0, T);

func real YConvexRightParabola(T);
real T;
return DeCasteljau(1.0, 2.0, 3.0, T);

func real ImplicitEquation(X, Y);
real X, Y;
return 5.0*X*X + Y - 5.0;

func real F(T);
real T;
real X, Y;
[X:= XConvexLeftParabola(T);
 Y:= YConvexRightParabola(T);
return ImplicitEquation(X, Y) + T;
];

func real Aitken(P0);
real P0;
real P1, P2, P1M0;
[P1:= F(P0);
 P2:= F(P1);
P1M0:= P1 - P0;
return P0 - P1M0 * P1M0 / (P2 - 2.0*P1 + P0);
];

func real SteffensenAitken(PInit, Tol, MaxIter);
real PInit, Tol; int MaxIter;
real P0, P;
int Iter;
[P0:= PInit;
P:= Aitken(P0);
Iter:= 1;
while abs(P-P0) > Tol and Iter < MaxIter do
        [P0:= P;
        P:= Aitken(P0);
        Iter:= Iter+1;
        ];
if abs(P-P0) > Tol then return NaN;
return P;
];

real T0, T, X, Y;
int  I;
[T0:= 0.0;
for I:= 0 to 10 do
        [Print("T0:= %1.1f : ", T0);
        T:= SteffensenAitken(T0, 0.00000001, 1000);
        if T = NaN then
                Print("No answer\n")
        else    [X:= XConvexLeftParabola(T);
                 Y:= YConvexRightParabola(T);
                if abs(ImplicitEquation(X, Y)) <= 0.000001 then
                        Print("Intersection at (%1.6f, %1.6f)\n", X, Y)
                else    Print("Spurious solution\n");
                ];
        T0:= T0 + 0.1;
        ];
]
Output:
T0:= 0.0 : No answer
T0:= 0.1 : Intersection at (0.881025, 1.118975)
T0:= 0.2 : No answer
T0:= 0.3 : No answer
T0:= 0.4 : No answer
T0:= 0.5 : No answer
T0:= 0.6 : No answer
T0:= 0.7 : No answer
T0:= 0.8 : No answer
T0:= 0.9 : Intersection at (-0.681025, 2.681025)
T0:= 1.0 : No answer