IMPLEMENTATION MODULE ApproxLib;
(*------------------------------------------------------------------------*)
(* Diese Modul stellt Routinen zur Ausgleichsbrechnung von Polynomen *)
(* und Splines zur Verfuegung *)
(* Module for fitting polynomals and splines to data sets and functions *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung *)
(* *)
(* 11.09.92, MRi: TschebKoeff eingefuegt *)
(* 24.07.95, MRi: Spline auf dynamische Speicherverwaltung umgeschrieben *)
(* 28.01.96, MRi: Vektorparameter der Routine FitPoly als VAR dekleriert *)
(* 05.05.96, MRi: In Ausgleichsgerade einige zus"atzliche Abfragen *)
(* eingef"ugt, so da3 diese Routine stabiler l"auft *)
(* 16.04.15, MRi: In Prozedur PolyKoeff TrVek auf Zeiger umgestellt, *)
(* Abfrage auf Laenge der Uebergabeparameter eingefuehrt *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 13.04.17, MRi: Erste Uebersetzung RatTsch1 von Algol60 nach Pascal *)
(* 14.04.17, MRi: Uebersetzung von RatTsch1 von Pascal nach Modula-2 *)
(* 20.02.18, MRi: In RatTsch1 bei der Berechnung von etabar (Label A5 - *)
(* A6) ein Korrekturterm (+ Small) eingefuehrt um *)
(* Divisionen durch Null zu verhindern *)
(* 21.02.18, MRi: Kosmetische Aenderungen an RatTsch1, Parametrisierung *)
(* 23.02.18, MRi: In Ausgleichsgerade SumY2 entfernt, RatTsch1 eingefuegt *)
(* 24.02.18, MRi: Modul von Ausgleich in ApproxLib umbenannt *)
(* 08.04.18, MRi: SplineKoeff und EvalSpline aus Fortran-Code (Ref. [1]) *)
(* erstellt. *)
(* 09.04.18, MRi: Prozedur Spline basiert nun auf SplineKoeff und *)
(* EvalSpline *)
(* 10.04.18, MRi: Prozedur EvalSplineDerivative hinzugefuegt *)
(*------------------------------------------------------------------------*)
(* Testroutine fue RatTsch1 in TstRatCheb (inklusieve GNUPLOT Ausgaben) *)
(* Testroutine fue Spline in TstSpline *)
(*------------------------------------------------------------------------*)
(* [1] Forsythe, George; Malcolm, Michael; Moler, Cleve; "Computer *)
(* Methods for Mathematical Computations", Prentice-Hall, *)
(* Upper Saddle River, US (1977) *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: ApproxLib.mod,v 1.6 2018/04/09 20:48:50 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM LowLong IMPORT intpart;
IMPORT Errors;
FROM Deklera IMPORT MaxDim;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
FROM LongMath IMPORT pi,sqrt,sin,cos;
IMPORT LMathLib;
FROM LMathLib IMPORT Funktion1;
FROM SpezFunkt1 IMPORT Clenshaw,Horner,RatPoly;
FROM MatLib IMPORT SumVek,ENorm;
FROM LinLib IMPORT HhLin;
IMPORT TIO;
CONST debug = FALSE;
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
PROCEDURE Ausgleichsgerade(VAR A,B : LONGREAL;
VAR dA,dB : LONGREAL;
VAR Xq,Yq : LONGREAL;
VAR R : LONGREAL;
dim : CARDINAL;
T : LONGREAL;
VAR X,Y : ARRAY OF LONGREAL);
CONST blas = TRUE;
eps = 2.6E-15; (* Zahl, unterhalb derer eine Summe *)
(* als Null betrachtet wird. *)
MaxXY = 8191;
VAR i : CARDINAL;
SumX,SumY,SumX2,SumXY : LONGREAL;
S,Sx,Sy,S1,S2,S3,S4,Dim : LONGREAL;
Nenner : LONGREAL;
XY : POINTER TO
ARRAY [0..MaxXY] OF LONGREAL;
dimflag : BOOLEAN;
(* Obergrenze in XY rein formal, vorsicht wenn Index"uber- *)
(* pr"ufung des "Ubersetzer eingeschaltet. Siehe Abfrage *)
(* auf dimflag im Quellcode. *)
BEGIN
Fehler:=FALSE;
IF (dim <= 2) THEN
ErrOut("Zuwenig Datenpunkte (Ausgleichsgerade)");
Fehler:=TRUE;
RETURN;
END;
dimflag := (dim > (MaxXY+1)); (* Abfrage wegen Index"uberschreitung *)
Dim:=VAL(LONGREAL,dim);
IF NOT blas OR dimflag THEN (* Scheller, aber weniger genau ! *)
SumX :=0.0; SumY :=0.0;
SumX2:=0.0; SumXY:=0.0;
FOR i:=0 TO dim-1 DO
SumX := SumX + X[i];
SumY := SumY + Y[i];
SumX2 := SumX2 + X[i]*X[i];
SumXY := SumXY + X[i]*Y[i];
END;
ELSE
SumX := SumVek(X,1,dim);
SumY := SumVek(Y,1,dim);
SumX2 := ENorm(X,1,dim); SumX2:=SumX2*SumX2;
ALLOCATE(XY,dim*TSIZE(LONGREAL));
FOR i:=0 TO dim-1 DO (* Das ist auch noch nicht so doll ! *)
XY^[i] := X[i]*Y[i];
END;
SumXY := SumVek(XY^,1,dim);
DEALLOCATE(XY,dim*TSIZE(LONGREAL));
END;
Xq := SumX / Dim;
Yq := SumY / Dim;
Nenner := SumX2 - (SumX*SumX / Dim);
IF (ABS(Nenner) < (eps * SumXY)) OR (Nenner = 0.0) THEN
(* Nenner Null, wenn alle X_i = 1 oder 0 *)
(* Dann sind die Eingabedaten aber nicht sinnvoll. *)
A:=MAX(LONGREAL); dA:=MAX(LONGREAL);
B:=MAX(LONGREAL); dB:=MAX(LONGREAL); R:=MAX(LONGREAL);
Fehlerflag:=
"Operation unbestimmt, Berechnung nicht moeglich (Ausgleichsgerade) !";
Fehler:=TRUE; ErrOut(Fehlerflag);
RETURN;
END;
A := (SumXY - Xq*SumY) / Nenner;
B := Yq - A*Xq;
S1:=0.0; S2:=0.0; S3:=0.0; S4:=0.0;
FOR i:=0 TO dim-1 DO
S1 := S1 + sqr(X[i] - Xq);
S2 := S2 + sqr(Y[i] - Yq);
S3 := S3 + (X[i] - Xq)*(Y[i] - Yq);
S4 := S4 + sqr(A*(X[i]) + B - Y[i]);
END;
Sx := sqrt(S1);
Sy := sqrt(S2);
S := sqrt(S4 / Dim);
IF (Sx # 0.0) THEN
dA := T*S / (sqrt((Dim-2.0) / Dim)*Sx);
ELSE
dA := 0.0;
END;
dB := T*S / sqrt (Dim-2.0);
IF (Sx*Sy # 0.0) THEN
R := S3 / (Sx*Sy);
ELSE
R := 0.0;
END;
END Ausgleichsgerade;
PROCEDURE TschebKoeff( U,O : LONGREAL;
VAR TKoeff : ARRAY OF LONGREAL;
Grad : CARDINAL;
Funkt : Funktion1);
(*----------------------------------------------------------------*)
(* C_i = {2 \oper N} \sum_{k=1}^N *)
(* f \left[ *)
(* \cos({\pi (k - 1/2) \over N} *)
(* \cos({\pi (j - 1) (k - 1/2) \over N} *)
(* \right] *)
(*----------------------------------------------------------------*)
VAR j,k : CARDINAL;
y,N,OmU,OpU,Sum,Fac,Jm1 : LONGREAL;
F : POINTER TO ARRAY [1..8192] OF LONGREAL;
BEGIN
Fehler:=TRUE;
IF (Grad+1 > HIGH(TKoeff)) THEN
Fehlerflag:="Koeffizientenvektor zu klein (ApproxLib.TschebKoeff) !";
Fehler:=TRUE; ErrOut(Fehlerflag);
RETURN;
END;
ALLOCATE(F,(Grad+1)*TSIZE(LONGREAL));
IF (F = NIL) THEN
Fehlerflag:="Kein Freispeicher vorhandden (ApproxLib.TschebKoeff) !";
Fehler:=TRUE; ErrOut(Fehlerflag);
RETURN;
END;
N := VAL(LONGREAL,Grad+1);
OmU := 0.5*(O - U);
OpU := 0.5*(O + U);
FOR k:=1 TO Grad+1 DO
y := cos(pi*(VAL(LONGREAL,k) - 0.5) / N);
F^[k] := Funkt(y*OmU + OpU);
END;
Fac := 2.0 / N;
FOR j:=1 TO Grad+1 DO
Sum:=0.0; Jm1 := VAL(LONGREAL,j-1);
FOR k:=1 TO Grad+1 DO
Sum:=Sum + F^[k]*cos(pi*Jm1 * ((VAL(LONGREAL,k) - 0.5) / N));
END;
TKoeff[j-1]:=Fac*Sum;
END;
DEALLOCATE(F,(Grad+1)*TSIZE(LONGREAL));
END TschebKoeff;
PROCEDURE PolyKoeff(VAR PKoeff : ARRAY OF LONGREAL; (* Polynomkoeffizienten *)
VAR TKoeff : ARRAY OF LONGREAL; (* Tschebytscheffkoeff. *)
U,O : LONGREAL; (* Intervall f"ur Polynom *)
Grad : CARDINAL);
VAR N,i,j,k : CARDINAL;
TrVek : POINTER TO ARRAY [0..8192] OF LONGREAL;
Fac,Konst,Zw : LONGREAL;
BEGIN
N := Grad + 1; (* Anzahl der Koeffizienten *)
IF (N > HIGH(TKoeff)) OR (N > HIGH(PKoeff)) THEN
Fehlerflag:="Koeffizientenvektoren zu klein (ApproxLib.PolyKoeff) !";
Fehler:=TRUE; ErrOut(Fehlerflag);
RETURN;
END;
IF (N > 8191) THEN
FatalError("Grad des Polynoms zu gro3 (ApproxLib.PolyKoeff) !");
END;
ALLOCATE(TrVek,N*TSIZE(LONGREAL));
FOR i:=0 TO Grad DO PKoeff[i]:=0.0; TrVek^[i]:=0.0; END;
PKoeff[0]:=TKoeff[Grad];
FOR j:=N-1 TO 2 BY -1 DO
FOR k:=N-j+1 TO 2 BY -1 DO
Zw:=PKoeff[k-1];
PKoeff[k-1]:= 2.0*PKoeff[k-2] - TrVek^[k-1];
TrVek^[k-1]:=Zw;
END;
Zw:=PKoeff[0];
PKoeff[0] := - TrVek^[0] + TKoeff[j-1];
TrVek^[0]:=Zw;
END; (* FOR j *)
FOR j:=N TO 2 BY -1 DO PKoeff[j-1] := PKoeff[j-2] - TrVek^[j-1]; END;
PKoeff[0] := - TrVek^[0] + 0.5*TKoeff[0];
Konst := 2.0 / (O - U); (* Transfomiere auf Intervall [U,O] *)
Fac := Konst;
FOR j:=2 TO N DO
PKoeff[j-1] := PKoeff[j-1] * Fac;
Fac := Fac * Konst;
END;
Konst := 0.5*(U + O);
FOR j:=1 TO N-1 DO
FOR k:=N-1 TO j BY -1 DO
PKoeff[k-1] := PKoeff[k-1] - Konst*PKoeff[k];
END;
END;
DEALLOCATE(TrVek,N*TSIZE(LONGREAL));
END PolyKoeff;
PROCEDURE RatKoeff( Funkt : Funktion1;
U,O : LONGREAL;
VAR Ak,Bk : ARRAY OF LONGREAL;
n,m : CARDINAL);
(*----------------------------------------------------------------*)
(* Nachsehen ob man nicht einen eren Algorithmus finden kann, *)
(* z.B. (durchsehen ob brauchbar) *)
(* *)
(* [1] Barrodale, I.; Mason, J.C.; "Two simple algorithms for *)
(* discrete rational approximation", *)
(* Math. Comp. 24, pp 877-891(1970) *)
(* [1] Van Barel, Marc; Bultheel, Adhemar; "A new approach to the *)
(* rational interpolation problem", *)
(* J. Comput. Appl. Math. 32, pp 281-289 (1990) *)
(* [2] Schock, Eberhard; "Pointwise Rational Approximation and *)
(* Iterative Methods for Ill-Posed Problems", *)
(* Num. Math. 54, pp 91 (1999) *)
(* [1] Van Barel, Marc; Bultheel, Adhemar; "A parallel algorithm *)
(* for discrete least squares rational approximation, *)
(* Num. Math. 63, pp. 99 (2008) *)
(*----------------------------------------------------------------*)
CONST NDim = MaxDim DIV 2;
VAR i : CARDINAL;
x,Delta,StdAbw : LONGREAL;
X,Y : ARRAY [0..NDim-1] OF LONGREAL;
BEGIN
x:=U; Delta := (O - U) / VAL(LONGREAL,NDim-2);
FOR i:=0 TO NDim-1 DO X[i]:=x; Y[i]:=Funkt(x); x:=x + Delta; END;
FitRat(X,Y,NDim,Ak,Bk,n,m,StdAbw);
END RatKoeff;
PROCEDURE FitTscheb(dim : CARDINAL;
Grad : CARDINAL;
VAR X,Y : ARRAY OF LONGREAL;
VAR TKoeff : ARRAY OF LONGREAL;
VAR U,O : LONGREAL;
VAR StdAbw : LONGREAL;
VAR IFehl : CARDINAL);
(*----------------------------------------------------------------*)
(* Die Tschebytscheffkoeffizienten werden als L"osung des *)
(* linearen Gleichungssystems F X - R = C erhalten *)
(* *)
(* F_{i,0} = 1, *)
(* F_{i,1} = x_j, *)
(* F_{i,j} = 2 F_{i,j-1} X_j - T{i,j-2} X_j \forall j > 1 *)
(*----------------------------------------------------------------*)
VAR i,j,ij : CARDINAL;
Sum,Ri,K1,K2 : LONGREAL;
F,C,X11 : POINTER TO ARRAY [0..8191] OF LONGREAL;
BEGIN
IFehl:=0; Fehler:=FALSE;
IF (dim < 2) THEN IFehl:=1; ELSIF (dim > 8191) THEN IFehl:=2; END;
IF (Grad < 1) THEN IFehl:=3; ELSIF (Grad > dim) THEN IFehl:=4; END;
IF ((Grad+1)*dim > 8191) THEN IFehl:=5; END;
CASE IFehl OF
1 : Fehlerflag:='Anzahl der Datenpunkte kleiner zwei.';|
2 : Fehlerflag:='Anzahl der Datenpunkte gr"o\3er als 8191.';|
3 : Fehlerflag:='Grad des T-Polynoms kleiner eins.';|
4 : Fehlerflag:='Grad des T-Polynoms gr"o\3er als Datenpunktanzahl.';|
5 : Fehlerflag:='Fehlergleichungssystem zu gro\3.';
ELSE
Fehlerflag:='Kein Fehler';
END;
IF (IFehl # 0) THEN Fehler:=TRUE; ErrOut(Fehlerflag); RETURN; END;
ALLOCATE(F,((Grad+1)*dim)*TSIZE(LONGREAL));
ALLOCATE(C,dim*TSIZE(LONGREAL));
ALLOCATE(X11,dim*TSIZE(LONGREAL));
IF (F = NIL) OR (C = NIL) OR (X11 = NIL) THEN
FatalError('Kein Freispeicher vorhanden (TschebFit) !');
END;
U:=X[0]; O:=X[0]; (* Transformiere X auf Intervall [-1,1] *)
FOR i:=1 TO dim-1 DO
IF (U > X[i]) THEN U := X[i]; END;
IF (O < X[i]) THEN O := X[i]; END;
END;
K1 := 0.5*(O + U);
K2 := 0.5*(O - U);
FOR i:=0 TO dim-1 DO X11^[i] := (X[i] - K1) / K2; END;
ij:=0;
FOR i:=1 TO dim DO (* Berechne die Fehlergleichungen *)
F^[ij] := 1.0; INC(ij);
IF (Grad > 0) THEN F^[ij] := X11^[i-1]; INC(ij); END;
FOR j:=2 TO Grad DO
F^[ij]:= 2.0*X11^[i-1] * F^[ij-1] - F^[ij-2]; INC(ij);
END;
END;
FOR i:=0 TO dim-1 DO C^[i]:= Y[i]; END;
HhLin(F^,TKoeff,C^,dim,Grad+1);
TKoeff[0]:=2.0*TKoeff[0];
IF Fehler THEN
Fehlerflag:='Keine Loesung gefunden (TschebKoeff) !';
IFehl:=6; StdAbw:=MAX(LONGREAL);
ELSE
Sum:=0.0; StdAbw:=0.0;
FOR i:=0 TO dim-1 DO
Ri:= Y[i] - Clenshaw(X11^[i],Grad,TKoeff,-1.0,1.0);
Sum:=Sum + Ri*Ri;
END;
IF (dim > Grad) THEN
StdAbw := sqrt(Sum / VAL(LONGREAL,dim - Grad));
END;
END;
DEALLOCATE(X11,dim*TSIZE(LONGREAL));
DEALLOCATE(C,dim*TSIZE(LONGREAL));
DEALLOCATE(F,((Grad+1)*dim)*TSIZE(LONGREAL));
END FitTscheb;
PROCEDURE FitPoly( dim : CARDINAL; (* Anzahl Datenpaare X,Y *)
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
Grad : CARDINAL; (* Grad des Polynoms *)
VAR PKoeff : ARRAY OF LONGREAL;(* Koeff. des Polynoms *)
VAR StdAbw : LONGREAL;
VAR IFehl : CARDINAL);
VAR i,j,ij : CARDINAL;
Sum,Ri,PotX : LONGREAL;
F,C : POINTER TO ARRAY [0..8192] OF LONGREAL;
BEGIN
IFehl:=0; Fehler:=FALSE;
IF (dim < 2) THEN IFehl:=1; ELSIF (dim > 8192) THEN IFehl:=2; END;
IF (Grad < 1) THEN IFehl:=3; ELSIF (Grad > dim) THEN IFehl:=4; END;
IF ((Grad+1)*dim > 8192) THEN IFehl:=5; END;
CASE IFehl OF
1 : Fehlerflag:='Anzahl der Datenpunkte kleiner zwei.';|
2 : Fehlerflag:='Anzahl der Datenpunkte gr"o3er als 8192.';|
3 : Fehlerflag:='Grad des Polynoms kleiner eins.';|
4 : Fehlerflag:='Grad des Polynoms gr"o3er als Datenpunktanzahl.';|
5 : Fehlerflag:='Fehlergleichungssystem zu gro3.';
ELSE
Fehlerflag:='Kein Fehler';
END;
IF (IFehl # 0) THEN Fehler:=TRUE; ErrOut(Fehlerflag); RETURN; END;
ALLOCATE(F,dim*(Grad+1)*TSIZE(LONGREAL));
ALLOCATE(C,dim*TSIZE(LONGREAL));
IF (F = NIL) OR (C = NIL) THEN
FatalError('Kein Freispeicher vorhanden (FitPoly) !');
END;
ij:=0; (* Berechne "gepacktes" Fehlergleichungssystem *)
FOR i:=0 TO dim-1 DO
C^[i]:=Y[i]; PotX:=1.0;
FOR j:=0 TO Grad DO F^[ij]:=PotX; PotX:=PotX * X[i]; INC(ij); END;
END;
HhLin(F^,PKoeff,C^,dim,Grad+1);
IF Fehler THEN
Fehlerflag:='Keine Loesung gefunden (FitPoly) !';
IFehl:=6; StdAbw:=MAX(LONGREAL);
ELSE
Sum:=0.0; StdAbw:=0.0;
FOR i:=0 TO dim-1 DO
Ri:= Y[i] - Horner(X[i],PKoeff,Grad);
Sum:=Sum + Ri*Ri;
END;
IF (dim > Grad) THEN
StdAbw := sqrt(Sum / VAL(LONGREAL,dim - Grad));
END;
END;
DEALLOCATE(C,dim*TSIZE(LONGREAL));
DEALLOCATE(F,dim*(Grad+1)*TSIZE(LONGREAL));
END FitPoly;
PROCEDURE FitRat(VAR X,Y : ARRAY OF LONGREAL;
dim : CARDINAL;
VAR Ak,Bk : ARRAY OF LONGREAL;
n,m : CARDINAL;
VAR StdAbw : LONGREAL);
VAR MaxK,Sum,d : LONGREAL;
Grad,i,j,ij,max : CARDINAL;
A,C,Koeff,PotX : POINTER TO ARRAY [0..8191] OF LONGREAL;
BEGIN
Grad := n + m + 2;
IF (Grad > dim+1) THEN
Fehler:=TRUE;
Fehlerflag:='Zuwenige Datenpunkte fuer die Interpolation (FitRat) !';
ErrOut(Fehlerflag);
RETURN;
END;
IF ((Grad)*dim > 8191) THEN
FatalError('Fehlergleichungssystem zu gro\3 (FitRat) !');
END;
IF (n > m) THEN max:=n; ELSE max:=m; END;
ALLOCATE(A,Grad*(dim+1)*TSIZE(LONGREAL));
ALLOCATE(C,(dim+1)*TSIZE(LONGREAL));
ALLOCATE(Koeff,Grad*TSIZE(LONGREAL));
ALLOCATE(PotX,(max+1)*TSIZE(LONGREAL));
IF (A = NIL) OR (C = NIL) OR (Koeff = NIL) OR (PotX = NIL) THEN
FatalError('Kein Freispeicher vorhanden (FitRat) !');
END;
ij:=0; PotX^[0]:=1.0;
FOR i:=0 TO dim-1 DO
FOR j:=1 TO max DO PotX^[j] := X[i]*PotX^[j-1]; END;
FOR j:=0 TO n DO A^[ij] := PotX^[j]; INC(ij); END;
FOR j:=0 TO m DO A^[ij] := - Y[i]*PotX^[j]; INC(ij); END;
C^[i]:=0.0;
END; (* FOR i *)
FOR i:=0 TO n+m+1 DO A^[ij]:=1.0; INC(ij); END;
C^[dim]:=1.0;
HhLin(A^,Koeff^,C^,(dim+1),Grad);
IF (Koeff^[n+1] # 0.0) THEN
MaxK := Koeff^[n+1];
ELSE
MaxK:=Koeff^[0];
FOR i:=1 TO Grad-1 DO
IF (ABS(Koeff^[i]) > ABS(MaxK)) THEN MaxK:=Koeff^[i]; END;
END;
END;
FOR i:=0 TO Grad-1 DO Koeff^[i] := Koeff^[i] / MaxK; END;
FOR i:=n TO HIGH(Ak) DO Ak[i] := 0.0; END;
FOR i:=m TO HIGH(Bk) DO Bk[i] := 0.0; END;
FOR i:=0 TO n DO Ak[i] := Koeff^[i]; END;
FOR i:=0 TO m DO Bk[i] := Koeff^[i+n+1]; END;
Sum:=0.0;
FOR i:=0 TO dim-1 DO
d := Y[i] - RatPoly(X[i],Ak,Bk,n,m);
Sum:=Sum + d*d;
END;
StdAbw:=0.0;
IF (dim > Grad) THEN
StdAbw := sqrt(Sum / VAL(LONGREAL,dim - Grad));
END;
DEALLOCATE(PotX,(max+1)*TSIZE(LONGREAL));
DEALLOCATE(Koeff,Grad*TSIZE(LONGREAL));
DEALLOCATE(C,(dim+1)*TSIZE(LONGREAL));
DEALLOCATE(A,Grad*(dim+1)*TSIZE(LONGREAL));
END FitRat;
PROCEDURE SplineKoeff(VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL);
VAR i : CARDINAL;
h : LONGREAL;
BEGIN
IF (N < 2) THEN RETURN END;
IF (N < 3) THEN
B[0]:=(Y[1] - Y[0]) / (X[1] - X[0]); (* Lineare Interpolation *)
C[0]:=0.0;
D[0]:=0.0;
B[1]:=B[0];
C[1]:=0.0;
D[1]:=0.0;
RETURN;
END;
(* preparation *)
D[0] := X[1] - X[0];
C[1] := (Y[1] - Y[0]) / D[0];
FOR i:=1 TO N-2 DO
D[i]:=X[i+1] - X[i];
B[i]:=2.0*(D[i-1] + D[i]);
C[i+1]:=(Y[i+1] - Y[i])/D[i];
C[i]:=C[i+1] - C[i];
END; (* FOR *)
(* end conditions *)
B[0] := - D[0];
B[N-1] := - D[N-2];
C[0] := 0.0;
C[N-1] := 0.0;
IF (N # 3) THEN
C[0] := C[2] / (X[3] - X[1] ) - C[1] / (X[2] - X[0] );
C[N-1] := C[N-2] / (X[N-1] - X[N-3]) - C[N-3] / (X[N-2] - X[N-4]);
C[0] := C[0 ]*sqr(D[0]) / (X[3] - X[0] );
C[N-1] := - C[N-1]*sqr(D[N-2]) / (X[N-1] - X[N-4]);
END; (* IF *)
(* forward elimination *)
FOR i:=1 TO N-1 DO
h := D[i-1] / B[i-1];
B[i] := B[i] - h*D[i-1];
C[i] := C[i] - h*C[i-1];
END;
(* back substitution *)
C[N-1]:=C[N-1] / B[N-1];
FOR i:=N-2 TO 0 BY -1 DO
C[i]:=(C[i] - D[i]*C[i+1]) / B[i];
END;
(* compute spline coefficients *)
B[N-1] := (Y[N-1] - Y[N-2])/D[N-2] + D[N-2]*(C[N-2] + 2.0*C[N-1]);
FOR i:=0 TO N-2 DO
B[i] := (Y[i+1] - Y[i])/D[i] - D[i]*(C[i+1] + 2.0*C[i]);
D[i] := (C[i+1] - C[i])/D[i];
C[i] := 3.0*C[i];
END;
C[N-1] := 3.0*C[N-1];
D[N-1] := D[N-2];
END SplineKoeff;
PROCEDURE EvalSpline( x : LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL) : LONGREAL;
VAR i,j,k : CARDINAL;
dx : LONGREAL;
BEGIN
(* if x is ouside the [X_0..X_{N-1}] interval take a boundary value *)
IF (x <= X[0]) THEN
RETURN Y[0];
END; (* IF *)
IF (x >= X[N-1]) THEN
RETURN Y[N-1];
END; (* IF *)
(* binary search for for i, such that X[i] <= x <= X[i+1] *)
i:=0;
j:=N;
WHILE (j > i+1) DO
k := (i + j) DIV 2;
IF (x < X[k]) THEN j:=k; ELSE i:=k; END;
END;
(* evaluate spline interpolation *)
dx := x - X[i];
RETURN Y[i] + dx*(B[i] + dx*(C[i] + dx*D[i]));
END EvalSpline;
PROCEDURE EvalSplineDerivative( x : LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL;
VAR y : LONGREAL;
VAR dy : LONGREAL;
VAR d2y : LONGREAL;
VAR iFehl : INTEGER);
VAR i,j,k : CARDINAL;
dx : LONGREAL;
BEGIN
iFehl := 0;
(* if x is ouside the [X_0..X_{N-1}] interval take a boundary value *)
dx := (X[1] - X[0]);
IF (x < (X[0]-0.25*dx)) THEN
y := Y[0];
dy := 0.0; d2y:= 0.0;
iFehl := 1;
RETURN;
END; (* IF *)
dx := (X[N-1] - X[N-2]);
IF (x > (X[N-1]+0.25*dx)) THEN
y := Y[N-1];
dy := 0.0; d2y:= 0.0;
iFehl := 1;
RETURN;
END; (* IF *)
(* binary search for for i, such that X[i] <= x <= X[i+1] *)
i:=0;
j:=N;
WHILE (j > i+1) DO
k := (i + j) DIV 2;
IF (x < X[k]) THEN j:=k; ELSE i:=k; END;
END;
(* evaluate interpolation including first and second derivative *)
dx := x - X[i];
y := Y[i] + dx*(B[i] + dx*(C[i] + D[i]*dx));
dy := B[i] + (2.0*C[i] + 3.0*D[i]*dx)*dx;
d2y := 2.0*C[i] + 6.0*D[i]*dx;
END EvalSplineDerivative;
PROCEDURE Spline(VAR Xa,Ya : ARRAY OF LONGREAL; (* Eingabevektoren *)
N : CARDINAL; (* Anzahl der Eingabewerte *)
X1,Xm : LONGREAL;
VAR Xn,Yn : ARRAY OF LONGREAL; (* Ausgabevektor *)
M : CARDINAL); (* Anzahl Splinepunkte *)
CONST MaxVek = MAX(INTEGER) DIV TSIZE(LONGREAL);
VAR i : CARDINAL;
DeltaX : LONGREAL;
B,C,D : POINTER TO ARRAY [0..MaxVek-1] OF LONGREAL;
BEGIN
IF (M-1 > HIGH(Xn)) OR (M-1 > HIGH(Yn)) THEN
FatalError("Felder Xn,Yn zu klein (ApproxLib.Spline) !");
END;
IF ((X1 >= Xa[0]) AND (Xm <= Xa[N-1])) THEN
ALLOCATE(B,N*TSIZE(LONGREAL));
ALLOCATE(C,N*TSIZE(LONGREAL));
ALLOCATE(D,N*TSIZE(LONGREAL));
IF (B = NIL) OR (C = NIL) OR (D = NIL) THEN
FatalError("Zuwenig Freispeicher vorhanden (ApproxLib.Spline) !");
END;
SplineKoeff(Xa,Ya,B^,C^,D^,N);
DeltaX := ABS(Xm - X1) / LFLOAT(M-1);
FOR i:=0 TO M-1 DO
Xn[i] := X1 + LFLOAT(i)*DeltaX;
Yn[i] := EvalSpline(Xn[i],Xa,Ya,B^,C^,D^,N);
END;
DEALLOCATE(D,N*TSIZE(LONGREAL));
DEALLOCATE(C,N*TSIZE(LONGREAL));
DEALLOCATE(B,N*TSIZE(LONGREAL));
ELSE
FatalError('Wert(e) au3erhalb des zulaessigen Bereichs (Spline).');
END;
END Spline;
PROCEDURE RatTsch1( nl,nr : CARDINAL;
ap,ep : LONGREAL;
eps : LONGREAL;
Fct : Funktion1;
W : Funktion1;
VAR A,B : ARRAY OF LONGREAL;
VAR etabar : LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Adjustments made to the original code *)
(* *)
(* For nsch please refer to the original article (bottom of page *)
(* 304, original value was 10, increase due to more computational *)
(* power of todays computers ;-) *)
(* For "genau" please refer to the original article page 305 in *)
(* the line before section "test results". The original value was *)
(* 1.0E-03 but had been decresed her due to the use of double *)
(* precision arithmetic. *)
(* Another deviation form the original is the use of the const *)
(* "Small" to avoid possible devisions by zero (see source code). *)
(* Replaced "1.5707963" by LMathLib.halfpi (see source code). *)
(*----------------------------------------------------------------*)
VAR l,r : INTEGER; (* Due to nl,nr beeing unsigned *)
CONST nsch = 20; (* Anzahl Mesh-Punkte *)
MaxIter = 10; (* max num. of iteration *)
genau = 10.0*LMathLib.MachEpsR4; (* Abbruchkriterium *)
Small = LMathLib.Small; (* correction term *)
PROCEDURE entier(x : LONGREAL) : INTEGER;
(*------------------------------------------------------*)
(* largest integer not greater than x, *)
(* x - 1 < entier(x) <= x *)
(*------------------------------------------------------*)
VAR i : INTEGER;
BEGIN
i := VAL(INTEGER,intpart(x));
IF (LFLOAT(i) > x) THEN DEC(i); END;
RETURN i;
END entier;
PROCEDURE sign(x : LONGREAL) : LONGREAL;
BEGIN
IF ((x > 0.0)) THEN
RETURN 1.0;
ELSIF ((x < 0.0)) THEN
RETURN -1.0;
ELSE
RETURN 0.0; (* x = 0 !!! *)
END;
END sign;
PROCEDURE Rat(VAR pp,qq : LONGREAL; (* uses A,B & l,r *)
z : LONGREAL);
VAR j : INTEGER;
p,q : LONGREAL;
BEGIN
p:=A[l];
FOR j:=l-1 TO 0 BY -1 DO p:=A[j] + p*z; END;
q:=B[r];
FOR j:=r-1 TO 0 BY -1 DO q:=B[j] + q*z; END;
IF (ABS(q) < ABS(p)*1.0E-06) THEN
(* Polerkennung / detection of poles ? *)
pp:=MAX(LONGREAL);
qq:=MAX(LONGREAL);
ELSE
pp:=p; qq:=q;
END;
END Rat;
PROCEDURE Epsn(z : LONGREAL) : LONGREAL;
VAR p,q : LONGREAL;
BEGIN
Rat(p,q,z);
IF (p = MAX(LONGREAL)) (* AND (q = MAX(LONGREAL)) *) THEN
Errors.ErrOut('*** error in Epsn ***');
RETURN MAX(LONGREAL);
ELSE
RETURN (Fct(z) - p/q)*W(z);
END;
END Epsn;
VAR i,j,k,n,m,ii,jj,mm : INTEGER;
it,ito,vz,wh : INTEGER;
alpha,min,max : LONGREAL;
dx,hs,q,ra : LONGREAL;
sigmaeta,sigmalbd : LONGREAL;
vh,z,zmax,zmin,zy : LONGREAL;
c,eta,f,g,h,qq,t,y : ARRAY [1..MaxL+MaxR+2] OF LONGREAL;
x : ARRAY [0..MaxL+MaxR+3] OF LONGREAL;
xa,xb : ARRAY [0..MaxL+MaxR] OF LONGREAL;
done,exit : BOOLEAN;
output : BOOLEAN; (* Print out some results ... *)
BEGIN
output := (iFehl < 0);
l := VAL(INTEGER,nl); r := VAL(INTEGER,nr);
ra:=0.0; (* Wg. Compilerwarnung *)
(* Computation of a first guess y[i] and x[i] for the zeros *)
(* and extreme points of the error function, respectively *)
(* A1: *)
n:=l + r;
ito:=2*n; IF (ito < MaxIter) THEN ito:=10; END;
it := 0; (* Iterationszaehler *)
x[1]:=ap; iFehl:= 0;
hs:=LMathLib.halfpi / LFLOAT(n+1); (* hs:=1.5707963 / FLOAT(n+1); *)
vh:=ep-ap; alpha:=0.0;
FOR i:=1 TO n+1 DO
q:= sin((LFLOAT(i) - 0.5)*hs);
ra := q*q*vh + ap;
y[i] := ra; (* :=:= *)
q:= sin(LFLOAT(i)*hs);
x[i+1]:=q*q*vh + ap;
q:= Fct(ra); f[i] := q; (* :=:= *)
q:= q * W(ra);
IF (ABS(q) > alpha) THEN alpha := ABS(q); END;
END;
x[n+2] := ep; alpha := 20.0*alpha*eps;
LOOP (* A2: *)
(* Calculation of a first trial rational approximation *)
(* by interpolation using the points (y[i],f[i]) *)
IF (l < r) THEN
FOR i:=1 TO n+1 DO
IF (f[i] = 0.0) THEN
iFehl := 2;
EXIT; (* GOTO neint *)
ELSE
f[i] := 1.0 / f[i];
END;
END;
i := l-1; j:=r
ELSE
i := r-1; j:=l ;
END;
FOR k:=1 TO n+1 DO
h[k] := f[k]; c[k] := -f[k]*y[k];
END;
FOR k:=1 TO n DO
FOR ii:=k+1 TO n+1 DO
max := h[ii] - h[k];
min := y[ii] - y[k];
hs := c[ii] - c[k];
IF (k <= j-i) THEN
h[ii] := max / min; c[ii] := hs / min;
ELSE
IF (ABS(hs) < (ABS(max) + ABS(min))*eps) THEN
iFehl := 2;
EXIT; (* GOTO neint *)
END;
h[ii] := max / hs; c[ii] := min / hs;
END;
END; (* ii *)
END; (* k *)
IF (0 <= i) THEN
xb[0] := h[n+1]; xa[0] := h[n] -c[n]*xb[0];
END;
FOR jj:=1 TO i DO
k := n+1-2*jj; xb[jj] := xb[jj-1];
xa[jj] := xa[jj-1] - c[k-1]*xb[jj];
FOR ii:=jj-1 TO 1 BY -1 DO
xb[ii] := xb[ii-1] - y[k ]*xb[ii] - c[k ]*xa[ii];
xa[ii] := xa[ii-1] - y[k-1]*xa[ii] - c[k-1]*xb[ii];
END;
xb[0] := h[k] - y[k ]*xb[0] - c[k ]*xa[0];
xa[0] := h[k-1] - y[k-1]*xa[0] - c[k-1]*xb[0];
END;
FOR k:=j-i TO 1 BY -1 DO
(* xa[j-k+1] := if k <= j then xa[j-k] else h[k]; *)
IF (k <= j) THEN xa[j-k+1] := xa[j-k]; ELSE xa[j-k+1] := h[k]; END;
FOR ii:=j-k TO 0 BY -1 DO
(*xa[ii] := (if ii = 0 then h[k] else xa[ii-1]) -
y[k]*xa[ii] - (if ii <= i then c[k]*xb[ii] else 0) *)
IF ((ii = 0)) THEN
IF ((ii <= i)) THEN
xa[ii] := h[k] - y[k]*xa[ii] - c[k]*xb[ii]
ELSE
xa[ii] := h[k] - y[k]*xa[ii];
END;
ELSE
IF ((ii <= i)) THEN
xa[ii] := xa[ii-1] - y[k]*xa[ii] - c[k]*xb[ii]
ELSE
xa[ii] := xa[ii-1] - y[k]*xa[ii];
END;
END;
END; (* ii *)
END; (* k *)
B[0] := 1.0;
IF (l < r) THEN
hs := 1.0 / xa[0]; A[0]:=hs; (* :=:= *)
FOR i:=1 TO r DO
B[i] := xa[i]*hs;
END;
FOR i:=1 TO l DO
A[i] := xb[i-1]*hs;
END;
ELSE
FOR i:=0 TO l DO
A[i] := xa[i];
END;
FOR i:=1 TO r DO
B[i] := xb[i-1];
END;
END;
(* Pointwise search of the extrema of the error function*)
(* A3: *)
x[0]:=ap; x[n+3] := ep; ii:=2;
vz := VAL(INTEGER,sign(Epsn(x[1])));
FOR i:=0 TO n+2 DO
vz := - vz;
IF (x[i] < x[i+1]) THEN
zmin := x[i+1];
z:=x[i]; zmax:=z; (* :=:= *)
max := Epsn(z )*LFLOAT(vz);
min := Epsn(zmin)*LFLOAT(vz);
dx:=(zmin - zmax) / LFLOAT(nsch);
mm:=0; m:=nsch;
FOR k:=0 TO nsch DO
hs := Epsn(z)*LFLOAT(vz);
IF (hs < min) THEN
min:=hs; m:=k; zmin:=z;
END;
IF (hs > max) THEN
max:=hs; mm:=k; zmax:=z;
END;
z:=z + dx;
END;
(* Computation of corrected extrema by quadratic interpolation *)
IF ((ap <= zmax-dx) AND (mm < nsch)) THEN
q := Epsn(zmax + dx);
zy := Epsn(zmax - dx);
ra:= (max*LFLOAT(vz*2) - q - zy)*2.0;
IF (ra # 0.0) THEN
zmax:=zmax + (q - zy)*dx/ra;
END;
END;
IF ((0 < m) AND (m < nsch)) THEN
q := Epsn(zmin + dx);
zy := Epsn(zmin - dx);
ra := (min*LFLOAT(vz*2) - q - zy)*2.0;
IF (ra # 0.0) THEN
zmin:=zmin + (q - zy)*dx/ra;
END;
END;
IF (i = 0) THEN
IF (zmax < zmin) THEN x[0]:=zmax; END;
x[1]:=zmin
ELSE
IF ((x[i-1] < zmax) AND (zmax < zmin)) THEN x[i] := zmax; END;
IF (x[i] < zmin) THEN x[i+1]:=zmin; END;
END;
IF (i = n+2) THEN
x[n+2] := zmax;
IF (zmax < zmin) THEN
IF (ABS(Epsn(zmax)) < ABS(Epsn(zmin))) THEN ii:=ii-1; END;
ELSIF (LFLOAT(vz)*Epsn(zmin) < LFLOAT(vz)*Epsn(x[n+1])) THEN
x[n+1] := zmin;
END;
END;
END;
END;
IF (ABS(Epsn(x[1])) < ABS(Epsn(x[0]))) THEN ii:=ii - 3;END;
IF (ii # 2) THEN
IF (ii = -2) THEN
jj := entier(LFLOAT(n+3) / 2.0); (* jj:=(n+3) DIV 2; *)
IF (n > 2) THEN
IF (ABS(Epsn(x[jj-1])) < ABS(Epsn(x[jj+1]))) THEN
jj:=jj - 1
ELSE
jj:=jj + 1;
END;
IF (ABS(Epsn(x[jj-1])) < ABS(Epsn(x[j+1]))) THEN
jj:=jj - 1;
END;
END;
FOR i:=jj TO 1 BY -1 DO
x[i] := x[i-1];
END;
FOR i:=jj+1 TO n+2 DO
x[i] := x[i+1];
END;
ELSE
i := entier(1.2 + LFLOAT(n+1)*LFLOAT(1-ii) / 2.0);
IF (ABS(Epsn(x[n + 3 - i + ii])) > ABS(Epsn(x[i]))) THEN
(* for k:=i step ii to i+ii*(n+1) do *)
(* x[k] := x[k+ii]; *)
k:=i;
WHILE (k < i+ii*(n+1)) DO
x[k] := x[k+ii];
k:=k + ii;
END;
END;
END;
END;
vz := VAL(INTEGER,sign(Epsn(x[1]))); wh := vz; (* :=:= *)
FOR i:=1 TO n+2 DO
vz:= -vz; zy := x[i];
Rat(ra,q,zy);
IF (ra = MAX(LONGREAL)) THEN
iFehl := iFehl + 2;
EXIT; (* GOTO pol *)
END;
qq[i] := q;
ra := ra / q; t[i] := ra; (* :=:= *)
q := LFLOAT(vz*wh) / W(zy); g[i] := q; (* :=:= *)
(* Korrektur: eta[i]:=(ra - Fct(zy)) / q; *)
eta[i]:=(ra - Fct(zy)) / (q + Small);
IF (i > 1) THEN
IF (sign(qq[i]) # sign(qq[i-1])) THEN ii := -4; END;
END;
END;
(* Computation of etabar *)
(* A4: *)
sigmaeta := 0.0; sigmalbd := 0.0; (* :=:= *)
FOR i:=1 TO n+2 DO
vh:= g[i]*qq[i]*qq[i];
FOR k:=1 TO n+2 DO
(* Korrektur: IF (k # i) THEN vh := vh / (x[i] - x[k]); END; *)
IF (k # i) THEN vh := vh / ((x[i] - x[k]) + Small); END;
END;
sigmaeta := vh * eta[i] + sigmaeta;
sigmalbd := vh + sigmalbd;
END;
(* Korrektur: etabar := sigmaeta / sigmalbd; *)
etabar := sigmaeta / (sigmalbd + Small);
(* Check of the accuracy obtained and of the number *)
(* of iterations. Preparations for the next iteration *)
(* A5: *)
(* FOR i:=1 TO n+2 DO *)
exit:=TRUE; done:=FALSE;
i:=1;
WHILE NOT done AND (i <= n+2) DO (* genau statt 1.0E-03 *)
IF (ABS(etabar - eta[i]) > ABS(etabar)*genau + alpha) THEN
it:=it + 1;
IF (it = ito) THEN
iFehl:=1;
IF (ii+3 > 0) THEN iFehl:=iFehl + 3; END; (* GOTO nenpruef *)
done:=TRUE;
ELSE
FOR k:=1 TO n+1 DO
j := k + entier(LFLOAT(2*k)/LFLOAT(n+2));
y[k] := x[j];
f[k] := t[j] + (etabar - eta[j])*g[j];
END;
exit:=FALSE;
done:=TRUE; (* GOTO A2; *)
END;
END;
INC(i);
END; (* WHILE *)
IF (exit) THEN EXIT END;
END; (* LOOP *)
(*
* nenpruef: IF ii+3 > 0 THEN GOTO nenul; END;
* GOTO A6;
* neint: iFehl := 2;
* nenul: iFehl := iFehl + 1;
* pol: iFehl := iFehl + 2;
* A6:
*)
IF (output) THEN
(* Output of the extrema x[i] and of the number *)
(* it of iterations needed can be done here *)
TIO.WrLn;
TIO.WrStr("RatTsch: Anzahl der Iterationen = "); TIO.WrInt(it,1);
TIO.WrStr("RatTsch: Fehlrcode = "); TIO.WrInt(iFehl,1);
TIO.WrLn;
END;
END RatTsch1;
END ApproxLib.