Steffensen's method: Difference between revisions
m
→{{header|Phix}}: use pygments
No edit summary |
m (→{{header|Phix}}: use pygments) |
||
| (43 intermediate revisions by 12 users not shown) | |||
Line 1:
{{task|Mathematics}}
[[wp:Steffensen's_method|Steffensen's method]] is a numerical method to find the roots of functions. It is similar to [[wp:Newton's_method|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 [[wp:Bézier_curve|Bézier curves]].
===Steffensen's method===
Line 43:
end
The function <code>steffensen_aitken</code> will find a fixed point for us, but how can we use that to find a root? In particular, suppose one wants to find <math>t</math> such that <math>g(t) = 0</math>. Then what one can do (and what we will try to do) is find a fixed point <math>p</math> of <math>f(t) = g(t) + t</math>.
===The problem we wish to solve===
Line 81:
fun y_convex_left_parabola (t : double) : double =
de_casteljau (1.0, 2.0, 3.0, t)
Plugging <math>x(t)</math> and <math>y(t)</math> into the implicit equation, and writing <math>f(t)</math> 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 <code>steffensen_aitken</code> 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 <math>{\it pinit}=t_{0}=0.0,0.1,0.2,\cdots,0.9,1.0</math>. Choose reasonable settings for tolerance and maximum iterations. (The [[#ATS|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 <math>t</math> into <math>x(t)</math> and <math>y(t)</math> to find the intersection point <math>(x,y)</math>. Check that that this answer is correct, by plugging <math>(x,y)</math> into the implicit equation. If it is not correct, print that Steffensen's method gave a spurious answer. Otherwise print the value of <math>(x,y)</math>.
(''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}}==
{{Trans|ATS}}
Uses a <code>UNION(VOID,REAL)</code> mode to represent the optional REAL value.
<syntaxhighlight lang="algol68">
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
</syntaxhighlight>
{{out}}
Same as the ATS sample:
<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
</pre>
=={{header|ATS}}==
<syntaxhighlight lang="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
</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
</pre>
=={{header|C}}==
{{trans|ATS}}
<syntaxhighlight lang="c">#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;
}</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|C++}}==
{{trans|C}}
<syntaxhighlight lang="cpp">#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;
}</syntaxhighlight>
{{out}}
<pre>
Same as C example.
</pre>
=={{header|EasyLang}}==
{{trans|C}}
<syntaxhighlight lang="easylang">
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
.
</syntaxhighlight>
=={{header|FreeBASIC}}==
{{trans|ATS}}
<syntaxhighlight lang="vb">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</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|Go}}==
{{trans|C}}
<syntaxhighlight lang="go">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
}
}</syntaxhighlight>
{{out}}
<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>
=={{header|Java}}==
{{trans|ATS}}
<syntaxhighlight lang="java">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;
}
}
}</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|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}}
<!--(phixonline)-->
<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}}==
{{trans|ATS}}
<syntaxhighlight lang="ratfor">
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
</syntaxhighlight>
{{out}}
Compiled with f2c (which accounts for the ink-saving numerals):
<pre> 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
</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>
=={{header|Wren}}==
{{trans|ATS}}
{{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="wren">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
}</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|XPL0}}==
{{trans|C}}
<syntaxhighlight lang "XPL0">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;
];
]</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>
| |||