IMPLEMENTATION MODULE Differ;
(*------------------------------------------------------------------------*)
(* Modul zur Berechnung von numerischen Gradienten *)
(* Module for calculation of numerical gradients *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 28.01.96, MRi: Ableit1 implementiert *)
(* 01.07.96, MRi: Durchsicht *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* Loeschen nicht genutzter Variablen *)
(* 03.10.16, MRi: Kosmetik, Ableit1p7 erstellt (extern) *)
(* 03.10.16, MRi: Prozedur Ableit1p7 Koeff gekuerzt und hier eingefuegt *)
(*------------------------------------------------------------------------*)
(* Testroutine in TestNumAbl *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: Differ.mod,v 1.4 2016/10/04 07:49:58 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE,ADR;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
FROM SpezFunkt1 IMPORT Horner;
FROM ApproxLib IMPORT FitPoly;
FROM LMathLib IMPORT Funktion1;
PROCEDURE Ableit1p5(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
VAR dY1 : ARRAY OF LONGREAL; (* <== Erste Ableitung *)
dH : LONGREAL; (* Schrittweite *)
n : CARDINAL); (* Anzahl Datenpunkte *)
VAR i : CARDINAL;
Koeff : LONGREAL;
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
BEGIN
Fehler := FALSE;
IF (n > 8192) THEN
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit1p5) !");
Fehler:=TRUE; RETURN;
END;
IF (n < 3) THEN
ErrOut("Zuwenig Datenpunkte (n < 3) (Differ.Ableit1p5) !");
Fehler:=TRUE; RETURN;
END;
IF (dH <= 0.0) THEN
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p5) !');
Fehler:=TRUE; RETURN;
END;
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
IF (Ytemp = NIL) THEN
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p5) !");
END;
Ytemp^[1] := (0.5 / dH) * (-3.0*Y[0] + 4.0*Y[1] - Y[2]);
IF (n > 3) THEN
Koeff := (1.0 / (12.0 * dH));
FOR i:=1 TO n-4 DO
Ytemp^[i+1] := Koeff*(- 3.0*Y[i-1] - 10.0*Y[i ] + 18.0*Y[i+1]
- 6.0*Y[i+2] + Y[i+3]);
END;
IF (n > 4) THEN
Ytemp^[n-2] := Koeff*(Y[n-5] - 8.0*Y[n-4] + 8.0*Y[n-2] - Y[n-1]);
ELSE
Ytemp^[2] := (0.5 / dH) * (-3.0*Y[1] + 4.0*Y[2] - Y[3]);
END;
END; (* IF n > 3 *)
Ytemp^[n-1] := (0.5 / dH) * (-Y[n-3] + Y[n-1]);
Ytemp^[n ] := (0.5 / dH) * (+Y[n-3] - 4.0*Y[n-2] + 3.0*Y[n-1]);
FOR i:=1 TO n DO dY1[i-1]:=Ytemp^[i]; END;
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
END Ableit1p5;
PROCEDURE Ableit1p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
VAR dY1 : ARRAY OF LONGREAL; (* <== Erste Ableitung *)
dH : LONGREAL; (* Schrittweite *)
N : CARDINAL); (* Anzahl Datenpunkte *)
VAR i : CARDINAL;
Koeff : LONGREAL;
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
BEGIN
Fehler := FALSE;
IF (N > 8192) THEN
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit1p7) !");
Fehler:=TRUE; RETURN;
END;
IF (N < 7) THEN
ErrOut("Zuwenig Datenpunkte (n < 3) (Differ.Ableit1p7) !");
Fehler:=TRUE; RETURN;
END;
IF (dH <= 0.0) THEN
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p7) !');
Fehler:=TRUE; RETURN;
END;
ALLOCATE(Ytemp,N*TSIZE(LONGREAL));
IF (Ytemp = NIL) THEN
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p7) !");
END;
(* Koeff. von http://web.media.mit.edu/~crtaylor/calculator.html *)
Ytemp^[1] := ( - 147.0*Y[0] + 360.0*Y[1] - 450.0*Y[2] + 400.0*Y[3]
- 225.0*Y[4] + 72.0*Y[5] - 10.0*Y[6]) / (60.0*dH);
Ytemp^[2] := ( - 10.0*Y[0] - 77.0*Y[1] + 150.0*Y[2] - 100.0*Y[3]
+ 50.0*Y[4] - 15.0*Y[5] + 2.0*Y[6]) / (60.0*dH);
Ytemp^[3] := ( 2.0*Y[0] - 24.0*Y[1] - 35.0*Y[2] + 80.0*Y[3]
- 30.0*Y[4] + 8.0*Y[5] - 1.0*Y[6]) / (60.0*dH);
Koeff := (1.0 / (60.0 * dH));
FOR i:=3 TO N-4 DO
Ytemp^[i+1] := Koeff*( - 1.0*Y[i-3] + 9.0*Y[i-2]
- 45.0*Y[i-1] + 45.0*Y[i+1]
- 9.0*Y[i+2] + 1.0*Y[i+3]);
END;
(* (* -3,-2,-1,0,1,2 *)
Ytemp^[N-2] := ( - 2.0*Y[N-6] + 15.0*Y[N-5] - 60.0*Y[N-4] + 20.0*Y[N-3]
+ 30.0*Y[N-2] - 3.0*Y[N-1]) / (60.0*dH); *)
Ytemp^[N-2] := ( 3.4203893866172544E+00*Y[N-7] -
2.7363115092938072E+01*Y[N-6] +
1.0261168159851840E+02*Y[N-5] -
2.7363115092937826E+02*Y[N-4] +
1.1971362853160657E+02*Y[N-3] +
8.2089345278808270E+01*Y[N-2] -
6.8407787732336830E+00*Y[N-1]
) / (2.0522336319702848E+02*dH);
(* (* -3,-2,-1,0,1 *)
Ytemp^[N-1] := ( - 1.0*Y[N-5] + 6.0*Y[N-4] - 18.0*Y[N-3] + 10.0*Y[N-2]
+ 3.0*Y[N-1]) / (12.0*dH); *)
Ytemp^[N-1] := (-2.9104222964425162E+00*Y[N-7] +
2.1828167223318920E+01*Y[N-6] -
7.2760557411062450E+01*Y[N-5] +
1.4552111482212514E+02*Y[N-4] -
2.1828167223318757E+02*Y[N-3] +
1.1205125841303632E+02*Y[N-2] +
1.4552111482212510E+01*Y[N-1]
) / (8.7312668893275070E+01*dH);
(* (* -3,-2,-1,0 *)
Ytemp^[N-0] := ( - 2.0*Y[N-4] + 9.0*Y[N-3] - 18.0*Y[N-2] + 11.0*Y[N-1])
/ (6.0*dH); *)
Ytemp^[N-0] := ( 1.8863848217683325E+00*Y[N-7] -
1.3581970716731817E+01*Y[N-6] +
4.2443658489786955E+01*Y[N-5] -
7.5455392870732334E+01*Y[N-4] +
8.4887316979573910E+01*Y[N-3] -
6.7909853583659130E+01*Y[N-2] +
2.7729856879994144E+01*Y[N-1]
) / (1.1318308930609855E+01*dH);
FOR i:=1 TO N DO dY1[i-1]:=Ytemp^[i]; END;
DEALLOCATE(Ytemp,N*TSIZE(LONGREAL));
END Ableit1p7;
PROCEDURE Ableit2p5(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
dH : LONGREAL; (* Schrittweite *)
n : CARDINAL); (* Anzahl Datenpunkte *)
VAR i : CARDINAL;
recH2 : LONGREAL; (* 1 / dH^2 *)
rec12H2 : LONGREAL; (* 1 / 12*dH^2 *)
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
BEGIN
Fehler := FALSE;
IF (n > 8192) THEN
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit2p5) !");
Fehler:=TRUE; RETURN;
END;
IF (n < 5) THEN
Fehler:=TRUE; Fehlerflag:='Zu wenige Datenpunkte (Ableit2p5)';
ErrOut(Fehlerflag); RETURN;
END;
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
IF (Ytemp = NIL) THEN
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p5) !");
END;
recH2 := 1.0 / (dH*dH);
Ytemp^[1] := recH2*(2.0*Y[0] - 5.0*Y[1] + 4.0*Y[2] - Y[3]);
Ytemp^[2] := recH2*(2.0*Y[1] - 5.0*Y[2] + 4.0*Y[3] - Y[4]);
rec12H2 := recH2 / 12.0;
FOR i:=2 TO n-3 DO
Ytemp^[i+1] := rec12H2*( - Y[i-2] + 16.0*Y[i-1] - 30.0*Y[i]
+ 16.0*Y[i+1] - Y[i+2]);
END;
Ytemp^[n-1]:= recH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
Ytemp^[n ]:= recH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
END Ableit2p5;
PROCEDURE Ableit2p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
dH : LONGREAL; (* Schrittweite *)
n : CARDINAL); (* Anzahl Datenpunkte *)
VAR i : CARDINAL;
recdH2 : LONGREAL; (* 1 / dH^2 *)
rec12dH2 : LONGREAL; (* 1 / 12*dH^2 *)
rec180dH2 : LONGREAL; (* 1 / 180*dH^2 *)
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
BEGIN
Fehler:=FALSE;
IF (n < 7) THEN
Fehler:=TRUE; Fehlerflag:='Zu wenige Datenpunkte (Ableit2p7)';
ErrOut(Fehlerflag); RETURN;
END;
IF (n > 8192) THEN
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit2p7) !");
Fehler:=TRUE; RETURN;
END;
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
IF (Ytemp = NIL) THEN
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p7) !");
END;
recdH2 := 1.0 / (dH*dH);
Ytemp^[1] := recdH2*(2.0*Y[0] - 5.0*Y[1] + 4.0*Y[2] - Y[3]);
Ytemp^[2] := recdH2*(2.0*Y[1] - 5.0*Y[2] + 4.0*Y[3] - Y[4]);
rec180dH2 := recdH2/180.0;
FOR i:=2 TO n-5 DO
Ytemp^[i+1]:= rec180dH2*(
- 13.0*Y[i-2] + 228.0*Y[i-1] - 420.0*Y[i]
+ 200.0*Y[i+1] + 15.0*Y[i+2] - 12.0*Y[i+3]
+ 2.0*Y[i+4]);
END;
rec12dH2 := recdH2 / 12.0;
Ytemp^[n-3] := rec12dH2*( - Y[n-6] + 16.0*Y[n-5]
- 30.0*Y[n-4] + 16.0*Y[n-3]
- Y[n-2]);
Ytemp^[n-2] := rec12dH2*( - Y[n-5] + 16.0*Y[n-4]
- 30.0*Y[n-3] + 16.0*Y[n-2]
- Y[n-1]);
Ytemp^[n-1]:= recdH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
Ytemp^[n ]:= recdH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
END Ableit2p7;
PROCEDURE AbleitN(VAR X,Y : ARRAY OF LONGREAL;
VAR dY : ARRAY OF LONGREAL;
ndata : CARDINAL;
Grad : CARDINAL;
NAbl : CARDINAL);
VAR i,k,o : CARDINAL;
dim : CARDINAL;
anfang,AltAnf : INTEGER;
PKoeff : POINTER TO ARRAY [0..8191] OF LONGREAL;
Xloc,Yloc : POINTER TO ARRAY [0..8191] OF LONGREAL;
dYloc : POINTER TO ARRAY [0..8191] OF LONGREAL;
StdAbw : LONGREAL;
IFehl : CARDINAL;
grad : CARDINAL;
BEGIN
(* n. Ableitung sonst nicht m"oglich. *)
IF (Grad < NAbl+1) THEN Grad:=NAbl+1; END;
dim := Grad + 2;
IF (ndata < dim) THEN
FatalError("Zuwenig Datenpunkte (Ableit1) !");
END;
ALLOCATE(dYloc ,ndata *TSIZE(LONGREAL));
ALLOCATE(PKoeff,(Grad+1)*TSIZE(LONGREAL));
IF (PKoeff = NIL) OR (dYloc = NIL) THEN
FatalError("Kein Freispeicher vorhanden (Ableit1) !");
END;
grad := Grad; (* wg. Compilerwarnung *)
AltAnf:=MAX(INTEGER);
FOR i:=0 TO ndata-1 DO
(* Bilde lokale Felder Xloc,Yloc passender L"ange *)
anfang := i - (dim DIV 2);
IF (anfang < 0) THEN
anfang:= 0;
ELSIF (anfang > VAL(INTEGER,(ndata - dim))) THEN
anfang:= ndata - dim;
END;
Xloc := ADR(X[anfang]);
Yloc := ADR(Y[anfang]);
IF (AltAnf # anfang) THEN
(* Fitte ein Polynom des Grades Grad an diese Subvektoren. *)
FitPoly(dim,Xloc^,Yloc^,Grad,PKoeff^,StdAbw,IFehl);
(* Berechne die n. Ableitung des Punkts (X_i,Y_i) aus den *)
(* Polynomkoeffizienten. *)
grad := Grad;
FOR o:=NAbl TO 1 BY -1 DO
FOR k:=1 TO grad DO
PKoeff^[k-1]:=VAL(LONGREAL,k)*PKoeff^[k];
END;
DEC(grad);
END;
END; (* IF *)
(* Auswertung der 1. Ableitun am Punkt X[i]. *)
dYloc^[i] := Horner(X[i],PKoeff^,grad);
AltAnf:=anfang;
END; (* FOR i *)
FOR i:=0 TO ndata-1 DO (* Kopiere L"osung *)
dY[i]:=dYloc^[i];
END;
DEALLOCATE(dYloc ,ndata *TSIZE(LONGREAL));
DEALLOCATE(PKoeff,(Grad+1)*TSIZE(LONGREAL));
END AbleitN;
PROCEDURE Ableit1und2( X : LONGREAL;
H : LONGREAL;
Funk : Funktion1;
VAR Y1,Y2 : LONGREAL);
VAR FmH,FpH : LONGREAL;
BEGIN
FpH := Funk(X+H);
FmH := Funk(X-H);
Y1 := (FpH - FmH) / (2.0*H);
Y2 := (FpH - 2.0*Funk(X) + FmH) / (H*H);
END Ableit1und2;
END Differ.