74 lines
1.5 KiB
Text
74 lines
1.5 KiB
Text
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;
|
|
];
|
|
]
|