98 lines
2.9 KiB
Plaintext
98 lines
2.9 KiB
Plaintext
#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
|