Steffensen's method: Difference between revisions
m
→{{header|Phix}}: use pygments
m (→{{header|Phix}}: use pygments) |
|||
| (13 intermediate revisions by 9 users not shown) | |||
Line 102:
(''The [[#ATS|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===
* [[Bézier_curves/Intersections|Bézier curves/Intersections]]
=={{header|ALGOL 68}}==
Line 490 ⟶ 493:
<syntaxhighlight lang="easylang">
s = 1 - t
c01 = s * c0 + t * c1
c12 = s * c1 + t * c2
.
.
.
.
return r +
.
p1m0 = p1 - p0
.
if abs (p - p0) < tol
.
p0 = p
.
return number "nan"
.
for i to 11
numfmt 1 0
write "t0 = " & t0 & " : "
numfmt 3 0
if t <> t
# nan
print "no answer"
else
if abs r <= 0.000001
print "intersection at (" & x & " " & y & ")"
Line 551:
.
</syntaxhighlight>
=={{header|FreeBASIC}}==
Line 718 ⟶ 717:
<pre>
Same as C example
</pre>
=={{header|Haxe}}==
{{Trans|ATS|via ALGOL 68}}
<syntaxhighlight lang="haxe">
// Steffensens's method - translated from the ATS sample via Algol 68
class OptionalReal
{
public var v:Float;
public var hasValue:Bool;
public function new( value:Float, notNull:Bool )
{
v = value;
hasValue = notNull;
};
}
class Main {
// Aitken's extrapolation
static function aitken( f:Float->Float // function Float to Float
, p0:Float // initial fixed point estimate
):Float
{
var p1 = f( p0 );
var p2 = f( p1 );
var p1m0 = p1 - p0;
return p0 - ( p1m0 * p1m0 ) / ( p2 - ( 2 * p1 ) + p0 );
};
// finds fixed point p such that f(p) = p
static function steffensenAitken( f:Float->Float // function Float to Float
, pinit:Float // initial estimate
, tol:Float // tolerance
, maxiter:Float // maximum number of iterations
):OptionalReal // return a Float, IF tolerance is met
{
var p0 = pinit;
var p = aitken( f, p0 ); // first iteration
var i = 1;
while( i < ( maxiter - 1 ) && Math.abs( p - p0 ) > tol ) {
p0 = p;
p = aitken( f, p0 );
i += 1;
}
if( Math.abs( p - p0 ) > tol ) {
return new OptionalReal( 0, false );
} else {
return new OptionalReal( p, true );
}
};
static function deCasteljau( c0:Float // control point coordinates (one axis)
, c1:Float
, c2:Float
, t :Float // the indep}ent parameter
) :Float // value of x(t) or y(t)
{
var s = 1 - t;
var c01 = ( s * c0 ) + ( t * c1 );
var c12 = ( s * c1 ) + ( t * c2 );
return ( s * c01 ) + ( t * c12 );
};
static function xConvexLeftParabola( t:Float ):Float { return deCasteljau( 2, -8, 2, t ); }
static function yConvexLeftParabola( t:Float ):Float { return deCasteljau( 1, 2, 3, t ); }
static function implicitEquation( x:Float, y:Float ):Float { return ( 5 * x * x ) + y - 5; }
static function f( t:Float ):Float // find fixed points of this function
{
var x = xConvexLeftParabola( t );
var y = yConvexLeftParabola( t );
return implicitEquation( x, y ) + t;
};
static function main()
{
var t0:Float = 0;
for( i in 1...12 ) {
Sys.print( 't0 = $t0 : ' );
var result = steffensenAitken( f, t0, 0.00000001, 1000 );
if( ! result.hasValue ) {
Sys.println( 'no answer' );
} else {
var x = xConvexLeftParabola( result.v );
var y = yConvexLeftParabola( result.v );
if( Math.abs( implicitEquation( x, y ) ) <= 0.000001 ) {
Sys.println( 'intersection at ($x, $y)' );
} else {
Sys.println( 'spurious solution' );
}
}
t0 += 0.1;
}
}
}
</syntaxhighlight>
{{out}}
<pre>
t0 = 0 : no answer
t0 = 0.1 : intersection at (0.881024967590665, 1.11897503240933)
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.68102496759051, 2.68102496759069)
t0 = 1 : no answer
</pre>
Line 805 ⟶ 917:
t0 = 1.0 : no answer
</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq'''
<syntaxhighlight lang="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
</syntaxhighlight>
{{output}}
<pre>
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
</pre>
=={{header|Julia}}==
{{trans|C}}
<syntaxhighlight lang="julia">""" 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()
</syntaxhighlight>{{out}}
<pre>
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
</pre>
=={{header|MATLAB}}==
{{trans|Julia}}
<syntaxhighlight lang="MATLAB}}">
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
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Nim}}==
{{trans|C}}
<syntaxhighlight lang="Nim">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
</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|Phix}}==
{{trans|C}}
<!--
<syntaxhighlight lang="phix">
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
</syntaxhighlight>
Output same as C
=={{header|Python}}==
{{trans|C}}
<syntaxhighlight lang="python">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()
</syntaxhighlight>{{out}} Output is the same as C.
=={{header|Raku}}==
{{trans|Go}}
<syntaxhighlight lang="raku" line># 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;
}</syntaxhighlight>
You may [https://ato.pxeger.com/run?1=jVTNbtswDL7slKdgDWe1m9izPAwommbFUOw2FMXOwwAllVat_skseXVQ5El26WF7mT1Cn2akFHluOwMzElkUyY_8SFo_fjb8pr2__9UamRw_vPit2xVwZW5EFYVyDuEmi-FuAgDlFoUcN0sIZeRE5gRrFC-gN2NlRhrUJ4SwIEUjTNtUJNpDNDmy6yt0z5M8teKMoBaT3WRCeWgjpBSVFtW7YUaqUgbfpi5wLXmnjGj-JhmFJJNdjDmw-ZCNTc6lRiCxTez2WhUC-EpjIoml8taCA6-uwILBaR_HhcEn3CDM40JZXIdJBmQ-m5G0G_D_R5yzM7jgF3BwgO6e-pU450i_-MrbKFxniL9mtOREfNARTRxTW1LzqMz6KMI_usIMVbhhsdu5Y-aP89iH7M7r6rvoPghpLnnDV3XBIxvLYw5zwn7NITmm1W7JcI-zdTgf1Zfr_wBiPcTrJziq3BRqrcz7by03qsYyd6jfDnDe0NB0-CMyW6wBHnhvGfV1egY0wnQ-mrorHRWKam6o01mKYy3rBj4z5sJsGlUZCKz2Do1SWZoTOJxmKZOHOziBYOG7ZtDk2XC_pCkyWAOEdg92nOHLjZSSqE2VplHxQ6j5FoKqxknVt6Jx-DsQhRa9if0i9oVbjvV4nLmfZoxOkzvWkxhOlz5v1sf2RZEQ4CoaLdbkBNxANEW2Uxl_qgIEsCA-1BMCnqbetI2qWw26LlqCCXqH_huj1swoEYZ3iLvS9jebv-H-AA Attempt This Online!]
=={{header|RATFOR}}==
Line 996 ⟶ 1,539:
t0 = .9 : intersection at ( -.681024968, 2.68102497)
t0 = 1. : no answer
</pre>
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
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
}
}
}
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
Line 1,002 ⟶ 1,624:
{{libheader|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.
<syntaxhighlight lang="
var aitken = Fn.new { |f, p0|
| |||