64 lines
1.4 KiB
Perl
64 lines
1.4 KiB
Perl
# 20240922 Perl programming solution
|
|
|
|
use strict;
|
|
use warnings;
|
|
sub aitken {
|
|
my ($f, $p0) = @_;
|
|
my $p2 = $f->( my $p1 = $f->($p0) );
|
|
return $p0 - ($p1 - $p0)**2 / ($p2 - 2 * $p1 + $p0);
|
|
}
|
|
|
|
sub steffensenAitken {
|
|
my ($f, $pinit, $tol, $maxiter) = @_;
|
|
my ($iter,$p) = (1, aitken($f, my $p0 = $pinit) );
|
|
|
|
while (abs($p - $p0) > $tol && $iter < $maxiter) {
|
|
$p = aitken($f, $p0 = $p);
|
|
$iter++
|
|
}
|
|
return abs($p - $p0) > $tol ? 'NaN' : $p;
|
|
}
|
|
|
|
sub deCasteljau {
|
|
my ($c0, $c1, $c2, $t) = @_;
|
|
my $s = 1.0 - $t;
|
|
return $s * ($s * $c0 + $t * $c1) + $t * ($s * $c1 + $t * $c2);
|
|
}
|
|
|
|
sub xConvexLeftParabola {
|
|
my ($t) = @_;
|
|
return deCasteljau(2.0, -8.0, 2.0, $t);
|
|
}
|
|
|
|
sub yConvexRightParabola {
|
|
my ($t) = @_;
|
|
return deCasteljau(1.0, 2.0, 3.0, $t);
|
|
}
|
|
|
|
sub implicitEquation {
|
|
my ($x, $y) = @_;
|
|
return 5.0 * ($x ** 2) + $y - 5.0;
|
|
}
|
|
|
|
sub f {
|
|
my ($t) = @_;
|
|
return implicitEquation(xConvexLeftParabola($t), yConvexRightParabola($t)) + $t;
|
|
}
|
|
|
|
for (my $t0 = 0.0, my $i = 0; $i < 11; ++$i) {
|
|
printf("t0 = %.1f : ", $t0);
|
|
my $t = steffensenAitken(\&f, $t0, 1e-8, 1000);
|
|
|
|
if ($t eq 'NaN') {
|
|
print "no answer\n";
|
|
} else {
|
|
my ($x, $y) = (xConvexLeftParabola($t), yConvexRightParabola($t));
|
|
if (abs(implicitEquation($x, $y)) <= 1e-6) {
|
|
printf("intersection at (%f, %f)\n", $x, $y);
|
|
} else {
|
|
print "spurious solution\n";
|
|
}
|
|
}
|
|
$t0 += 0.1;
|
|
}
|