IMPLEMENTATION MODULE SpezFunkt1;
(*------------------------------------------------------------------------*)
(* Stellt einige "Sonderfunktionen" fuer Polynome und Nullstellen- *)
(* berechnungen bereit. *)
(* Provides some special procedure for the evaluation of polynominals and *)
(* for the calculation of roots. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 17.07.93, MRi: Durchsicht *)
(* 07.10.15, MRi: Umbenennen von Funkt nach SpezFunkt1 *)
(* 02.01.16, MRi: Prozedur EvalCheb1Coeff eingefuegt. *)
(* 04.01.16, MRi: Prozedur ChebPoly1Sum eingefuegt *)
(* 29.06.16, MRi: Prozedur Nullstellen erweitert *)
(* 29.06.16, MRi: Prozedur Regula um Test auf Ueberlauf erweitert *)
(* 02.10.16, MRi: Prozedur QuadGl und KubGl aus ISO COMPLEX Typ umge- *)
(* stellt *)
(* 19.09.17, MRi: Korrektur in QuadGl fuer X1re = 0 *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SpezFunkt1.mod,v 1.6 2016/10/12 09:22:26 mriedl Exp mriedl $ *)
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut;
FROM LongMath IMPORT pi,sqrt,cos,arccos,exp,ln;
FROM LongComplexMath IMPORT zero,scalarMult,conj;
IMPORT LongComplexMath;
FROM LMathLib IMPORT MachEps,Small,sqrt3,
sign2,sqr,CardPot,MinCard,Funktion1;
PROCEDURE Horner(VAR x : LONGREAL;
VAR Koeff : ARRAY OF LONGREAL;
n : CARDINAL) : LONGREAL;
VAR y : LONGREAL;
i : CARDINAL;
BEGIN
y:=Koeff[n]; FOR i:=n-1 TO 0 BY -1 DO y:=y*x + Koeff[i]; END;
RETURN y;
END Horner;
PROCEDURE Clenshaw(x : LONGREAL; (* Funktionswert *)
n : CARDINAL; (* Grad d. T-Entwicklung *)
VAR C : ARRAY OF LONGREAL;
U,O : LONGREAL) : LONGREAL;
VAR i : CARDINAL;
y,y2,D2,D3,Zw : LONGREAL;
BEGIN
y := (2.0*x - O - U) / (O - U); y2 := 2.0*y;
D2:=0.0; D3:=0.0;
FOR i:=n TO 1 BY -1 DO Zw:=D2; D2:=y2*D2 - D3 + C[i]; D3:=Zw; END;
RETURN y*D2 - D3 + 0.5*C[0];
END Clenshaw;
PROCEDURE EvalCheb1Coeff( x : LONGREAL;
VAR T : ARRAY OF LONGREAL;
n : CARDINAL);
VAR i : CARDINAL;
BEGIN
T[0]:=1.0;
IF (n >= 1) THEN T[1]:=x; END;
IF (n >= 2) THEN
FOR i:=2 TO n DO
T[i]:=2.0*x*T[i-1] - T[i-2];
END;
END;
END EvalCheb1Coeff;
PROCEDURE ChebPoly1Sum( n : CARDINAL;
x : LONGREAL;
VAR A : ARRAY OF LONGREAL) : LONGREAL;
VAR k : CARDINAL;
h,r,s,tx : LONGREAL;
BEGIN
IF (n = 0) THEN
RETURN A[0];
ELSIF (n = 1) THEN
RETURN A[0] + A[1]*x;
ELSE
tx := x + x;
r := A[n];
h := A[n-1] + r*tx;
FOR k:=n-2 TO 1 BY -1 DO
s := r;
r := h;
h := A[k] + r*tx - s;
END;
RETURN A[0] - r + h*x;
END;
END ChebPoly1Sum;
PROCEDURE RatPoly( x : LONGREAL;
VAR A,B : ARRAY OF LONGREAL;
n,m : CARDINAL) : LONGREAL;
VAR i : CARDINAL;
Zaehler,Nenner : LONGREAL;
BEGIN
Zaehler:=A[n];
FOR i:=n-1 TO 0 BY -1 DO Zaehler := x*Zaehler + A[i]; END;
Nenner :=B[m];
FOR i:=m-1 TO 0 BY -1 DO Nenner := x*Nenner + B[i]; END;
IF (ABS(Nenner) > 1.0E-34) THEN
RETURN Zaehler / Nenner;
ELSE
RETURN MAX(LONGREAL);
END;
END RatPoly;
PROCEDURE QuadGl( A,B,C : LONGREAL;
VAR X1,X2 : LONGCOMPLEX);
VAR Nenner,Dis : LONGREAL;
X1re,X1im,X2re,X2im : LONGREAL;
BEGIN
Dis := B*B - 4.0*A*C;
Nenner := 2.0*A;
X1re := - B / Nenner;
X2re := X1re;
IF (Dis < 0.0) THEN (* Komplexe Wurzeln *)
X1im := - sqrt(-Dis) / Nenner;
X2im := - X1im;
ELSE (* Reelle Wurzeln *)
X1im := 0.0;
X2im := 0.0;
IF (B < 0.0) THEN
X1re := X1re + sqrt(Dis) / Nenner;
ELSE
X1re := X1re - sqrt(Dis) / Nenner;
END;
IF (X1re = 0.0) THEN
X2re:=0.0; (* ??? *)
ELSE
X2re := C / (A*X1re);
END;
END;
X1 := CMPLX(X1re,X1im);
X2 := CMPLX(X2re,X2im);
END QuadGl;
PROCEDURE KubGl( A,B,C,D : LONGREAL;
VAR X1,X2,X3 : LONGCOMPLEX);
PROCEDURE Root3(x : LONGREAL) : LONGREAL;
VAR sign : LONGREAL;
BEGIN
sign:=1.0; IF (x < 0.0) THEN sign:=-1.0; END;
RETURN sign*exp(ln(ABS(x))/3.0);
END Root3;
CONST sqrt3 = 1.7320508075688772;
W1 = CMPLX(-0.5, 0.5*sqrt3);
W2 = CMPLX(-0.5,-0.5*sqrt3);
VAR P,Q,Z,T,Phi : LONGREAL;
Z1,Z2,Cp,Px,A0,A3 : LONGREAL;
V,U : LONGCOMPLEX;
Y1,Y2,Y3 : LONGCOMPLEX;
BEGIN
A0 := A; A := B/A0; B := C/A0; C := D/A0;
P := - sqr(A)/9.0 + B/3.0;
Q := A*(sqr(A)/27.0 - B/6.0) + C/2.0;
Z := sqr(Q) + CardPot(P,3);
A3 := A / 3.0;
IF (ABS(Z) < 1.0E-15) THEN
T := Root3(-Q);
Y1 := CMPLX(2.0*T - 0.0, 0.0);
Y2 := CMPLX( - T - 0.0, 0.0);
Y3 := CMPLX( - T - 0.0, 0.0);
ELSIF (Z > 0.0) THEN
Z1 := - Q + sqrt(Z);
Z2 := - Q - sqrt(Z);
U := CMPLX(Root3(Z1),0.0);
V := CMPLX(Root3(Z2),0.0);
Y1 := U + V;
Y2 := W1*U + W2*V;
Y3 := W2*U + W1*V;
ELSIF (Z < 0.0) THEN
Cp := - Q / sqrt(-CardPot(P,3));
Phi := arccos(Cp);
Px :=-2.0*sqrt(-P);
Y1 := CMPLX(-Px*cos( Phi /3.0), 0.0);
Y2 := CMPLX( Px*cos((Phi + pi)/3.0), 0.0);
Y3 := CMPLX( Px*cos((Phi - pi)/3.0), 0.0);
ELSE (* Unfug wg. Compilerwarnung *)
Y1 := zero; Y2 := zero; Y3 := zero;
END;
X1 := Y1 - CMPLX(A3,0.0);
X2 := Y2 - CMPLX(A3,0.0);
X3 := Y3 - CMPLX(A3,0.0);
END KubGl;
PROCEDURE CubicEq( A,B,C,D : LONGCOMPLEX; (* coeffs of the polynomial *)
VAR X1,X2,X3 : LONGCOMPLEX);
CONST Zeta1 = CMPLX(-0.5, 0.5*sqrt3);
Zeta2 = CMPLX(-0.5, -0.5*sqrt3);
third = 1.0 / 3.0;
VAR P,Q,P2,A2,A0,s1,s2,delta : LONGCOMPLEX;
BEGIN
A0 := A; A := - B/A0; B := C/A0; C := - D/A0;
A2 := A*A;
P := A*(scalarMult(2.0,A2) - scalarMult(9.0,B) + scalarMult(27.0,C));
Q := A2 - scalarMult(3.0,B);
P2 := P*P;
delta := LongComplexMath.sqrt(P2 - scalarMult(4.0,(Q*Q*Q)));
IF (RE(conj(P)*delta) >= 0.0) THEN
s1 := LongComplexMath.power(scalarMult(0.5,(P + delta)),third);
ELSE
s1 := LongComplexMath.power(scalarMult(0.5,(P - delta)),third);
END;
IF (s1 = zero) THEN
s2 := zero;
ELSE
s2 := Q / s1;
END;
X1 := scalarMult(third, A + s1 + s2);
X2 := scalarMult(third, A + s1*Zeta2 + s2*Zeta1);
X3 := scalarMult(third, A + s1*Zeta1 + s2*Zeta2);
END CubicEq;
PROCEDURE Regula(VAR Funkt : Funktion1; (* y:=Funkt(x), mit x,y : LONGREAL *)
Xu,Xo : LONGREAL; (* Startwerte fuer die Nullstelle *)
fXu,fXo : LONGREAL; (* Funktionswerte von Xu,Xo. *)
VAR X0 : LONGREAL; (* Nullstelle von Funkt *)
genau : LONGREAL); (* Geforderte Genauigkeit *)
CONST MaxIter = 32;
VAR Iter : CARDINAL;
Nenner : LONGREAL;
BEGIN
Iter:=0; Fehler:=FALSE;
LOOP
INC(Iter);
IF (Iter > MaxIter) THEN
Fehler:=TRUE; ErrOut('Iter > Maxiter (Regula).');
RETURN;
END;
Nenner := (fXo - fXu);
IF (ABS(Nenner) < Small) THEN
Nenner := sign2(Small,Nenner); (* resultiert in grosser Zahl *)
END;
X0 := Xo - fXo*(Xo - Xu) / Nenner;
IF ((ABS(X0 - Xo)) < genau*ABS(X0)) THEN RETURN END;
fXu:=fXo; fXo:=Funkt(X0);
Xu:=Xo; Xo:=X0;
END;
END Regula;
PROCEDURE NullStellen(u,o : LONGREAL; (* Intervallgrenzen *)
dH : LONGREAL; (* Schrittweite *)
f : Funktion1;
VAR X0 : ARRAY OF LONGREAL; (* Nullstellen von f *)
VAR nx0 : CARDINAL; (* Anzahl der Nullstellen *)
i : CARDINAL; (* Anzahl zu suchender Nullst. *)
genau : LONGREAL);
VAR xu,xo,fxu,fxo : LONGREAL;
max : CARDINAL;
BEGIN
IF (i = 0) THEN max:=HIGH(X0)+1 ELSE max:=MinCard(HIGH(X0)+1,i) END;
fxo:=f(u); xo:=u; nx0:=0;
IF (ABS(fxo) < genau) THEN
X0[0]:=u; INC(nx0);
IF (max = 1) THEN RETURN; END;
END;
REPEAT
xu:=xo; xo:=xo + dH;
fxu:=fxo; fxo:=f(xo);
IF (fxo*fxu < 0.0) THEN (* Vorzeichenwechsel *)
Regula(f,xu,xo,fxu,fxo,X0[nx0],genau);
INC(nx0);
ELSIF ((ABS(fxo)-MachEps) < genau) THEN (* Funktionswert ist Null *)
X0[nx0]:=xo;
INC(nx0);
xo:=xo + dH;
END;
IF (nx0 = max) THEN RETURN; END;
UNTIL (xo >= o);
Fehler := (nx0 = 0);
IF Fehler THEN
Fehlerflag:='Keine Nullstellen im Intervall.';
END;
END NullStellen;
END SpezFunkt1.