Data update

Task/Steffensens-method/MATLAB/steffensens-method.m

@@ -0,0 +1,83 @@
+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

Task/Steffensens-method/Scala/steffensens-method.scala

@@ -0,0 +1,59 @@
+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
+    }
+  }
+}