How to resolve the algorithm Steffensen's method step by step in the Raku programming language

Published on 12 May 2024 09:40 PM
#Raku

How to resolve the algorithm Steffensen's method step by step in the Raku programming language

Table of Contents

Problem Statement

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.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Steffensen's method step by step in the Raku programming language

Source code in the raku programming language

# 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 also check:How to resolve the algorithm 100 prisoners step by step in the Commodore BASIC programming language
You may also check:How to resolve the algorithm Last Friday of each month step by step in the V (Vlang) programming language
You may also check:How to resolve the algorithm Dijkstra's algorithm step by step in the Ring programming language
You may also check:How to resolve the algorithm Password generator step by step in the Haskell programming language
You may also check:How to resolve the algorithm Balanced brackets step by step in the Déjà Vu programming language