|
a/Differ.mod |
|
b/Differ.mod |
|
... |
|
... |
10 |
(* 01.07.96, MRi: Durchsicht *)
|
10 |
(* 01.07.96, MRi: Durchsicht *)
|
11 |
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
|
11 |
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
|
12 |
(* Loeschen nicht genutzter Variablen *)
|
12 |
(* Loeschen nicht genutzter Variablen *)
|
13 |
(* 03.10.16, MRi: Kosmetik, Ableit1p7 erstellt (extern) *)
|
13 |
(* 03.10.16, MRi: Kosmetik, Ableit1p7 erstellt (extern) *)
|
14 |
(* 03.10.16, MRi: Prozedur Ableit1p7 Koeff gekuerzt und hier eingefuegt *)
|
14 |
(* 03.10.16, MRi: Prozedur Ableit1p7 Koeff gekuerzt und hier eingefuegt *)
|
|
|
15 |
(* 02.02.19, MRi: TSIZE durch SIZE ersetztt, Kosmetik in Ableit1p7 *)
|
15 |
(*------------------------------------------------------------------------*)
|
16 |
(*------------------------------------------------------------------------*)
|
16 |
(* Testroutine in TestNumAbl *)
|
17 |
(* Testroutine in TestNumAbl *)
|
17 |
(*------------------------------------------------------------------------*)
|
18 |
(*------------------------------------------------------------------------*)
|
18 |
(* Implementation : Michael Riedl *)
|
19 |
(* Implementation : Michael Riedl *)
|
19 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
20 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
20 |
(*------------------------------------------------------------------------*)
|
21 |
(*------------------------------------------------------------------------*)
|
21 |
|
22 |
|
22 |
(* $Id: Differ.mod,v 1.4 2016/10/04 07:49:58 mriedl Exp mriedl $ *)
|
23 |
(* $Id: Differ.mod,v 1.4 2016/10/04 07:49:58 mriedl Exp mriedl $ *)
|
23 |
|
24 |
|
24 |
FROM SYSTEM IMPORT TSIZE,ADR;
|
25 |
FROM SYSTEM IMPORT ADR;
|
25 |
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
|
26 |
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
|
26 |
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
|
27 |
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
|
27 |
FROM SpezFunkt1 IMPORT Horner;
|
28 |
FROM SpezFunkt1 IMPORT Horner;
|
28 |
FROM ApproxLib IMPORT FitPoly;
|
29 |
FROM ApproxLib IMPORT FitPoly;
|
29 |
FROM LMathLib IMPORT Funktion1;
|
30 |
FROM LMathLib IMPORT Funktion1;
|
|
... |
|
... |
49 |
END;
|
50 |
END;
|
50 |
IF (dH <= 0.0) THEN
|
51 |
IF (dH <= 0.0) THEN
|
51 |
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p5) !');
|
52 |
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p5) !');
|
52 |
Fehler:=TRUE; RETURN;
|
53 |
Fehler:=TRUE; RETURN;
|
53 |
END;
|
54 |
END;
|
54 |
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
55 |
ALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
55 |
IF (Ytemp = NIL) THEN
|
56 |
IF (Ytemp = NIL) THEN
|
56 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p5) !");
|
57 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p5) !");
|
57 |
END;
|
58 |
END;
|
58 |
Ytemp^[1] := (0.5 / dH) * (-3.0*Y[0] + 4.0*Y[1] - Y[2]);
|
59 |
Ytemp^[1] := (0.5 / dH) * (-3.0*Y[0] + 4.0*Y[1] - Y[2]);
|
59 |
IF (n > 3) THEN
|
60 |
IF (n > 3) THEN
|
|
... |
|
... |
67 |
ELSE
|
68 |
ELSE
|
68 |
Ytemp^[2] := (0.5 / dH) * (-3.0*Y[1] + 4.0*Y[2] - Y[3]);
|
69 |
Ytemp^[2] := (0.5 / dH) * (-3.0*Y[1] + 4.0*Y[2] - Y[3]);
|
69 |
END;
|
70 |
END;
|
70 |
END; (* IF n > 3 *)
|
71 |
END; (* IF n > 3 *)
|
71 |
Ytemp^[n-1] := (0.5 / dH) * (-Y[n-3] + Y[n-1]);
|
72 |
Ytemp^[n-1] := (0.5 / dH) * (-Y[n-3] + Y[n-1]);
|
|
|
73 |
Ytemp^[n ] := (0.5 / dH) * ( Y[n-3] - 4.0*Y[n-2] + 3.0*Y[n-1]);
|
|
|
74 |
(***************
|
|
|
75 |
Original Line :
|
72 |
Ytemp^[n ] := (0.5 / dH) * (+Y[n-3] - 4.0*Y[n-2] + 3.0*Y[n-1]);
|
76 |
Ytemp^[n ] := (0.5 / dH) * (+Y[n-3] - 4.0*Y[n-2] + 3.0*Y[n-1]);
|
|
|
77 |
|
|
|
78 |
HINT: Here the "+Y[n-3]" crashes the GNU M2 compiler ...
|
|
|
79 |
if I remove the "+" the crash is gone (even space around the
|
|
|
80 |
"+" does not seem to help ...
|
|
|
81 |
***************)
|
73 |
FOR i:=1 TO n DO dY1[i-1]:=Ytemp^[i]; END;
|
82 |
FOR i:=1 TO n DO dY1[i-1]:=Ytemp^[i]; END;
|
74 |
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
83 |
DEALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
75 |
END Ableit1p5;
|
84 |
END Ableit1p5;
|
76 |
|
85 |
|
77 |
PROCEDURE Ableit1p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
86 |
PROCEDURE Ableit1p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
78 |
VAR dY1 : ARRAY OF LONGREAL; (* <== Erste Ableitung *)
|
87 |
VAR dY1 : ARRAY OF LONGREAL; (* <== Erste Ableitung *)
|
79 |
dH : LONGREAL; (* Schrittweite *)
|
88 |
dH : LONGREAL; (* Schrittweite *)
|
80 |
N : CARDINAL); (* Anzahl Datenpunkte *)
|
89 |
N : CARDINAL); (* Anzahl Datenpunkte *)
|
81 |
|
90 |
|
82 |
VAR i : CARDINAL;
|
91 |
VAR i : CARDINAL;
|
83 |
Koeff : LONGREAL;
|
92 |
Koeff,dH60 : LONGREAL;
|
84 |
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
|
93 |
Ytemp : POINTER TO ARRAY [1..8192] OF LONGREAL;
|
85 |
BEGIN
|
94 |
BEGIN
|
86 |
Fehler := FALSE;
|
95 |
Fehler := FALSE;
|
87 |
IF (N > 8192) THEN
|
96 |
IF (N > 8192) THEN
|
88 |
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit1p7) !");
|
97 |
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit1p7) !");
|
89 |
Fehler:=TRUE; RETURN;
|
98 |
Fehler:=TRUE; RETURN;
|
|
... |
|
... |
94 |
END;
|
103 |
END;
|
95 |
IF (dH <= 0.0) THEN
|
104 |
IF (dH <= 0.0) THEN
|
96 |
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p7) !');
|
105 |
ErrOut('Unzulaessige Schrittweite (Differ.Ableit1p7) !');
|
97 |
Fehler:=TRUE; RETURN;
|
106 |
Fehler:=TRUE; RETURN;
|
98 |
END;
|
107 |
END;
|
99 |
ALLOCATE(Ytemp,N*TSIZE(LONGREAL));
|
108 |
ALLOCATE(Ytemp,N*SIZE(LONGREAL));
|
100 |
IF (Ytemp = NIL) THEN
|
109 |
IF (Ytemp = NIL) THEN
|
101 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p7) !");
|
110 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit1p7) !");
|
102 |
END;
|
111 |
END;
|
103 |
|
112 |
|
104 |
(* Koeff. von http://web.media.mit.edu/~crtaylor/calculator.html *)
|
113 |
(* Koeff. von http://web.media.mit.edu/~crtaylor/calculator.html *)
|
105 |
|
114 |
|
|
|
115 |
dH60 := 1.0 / (60.0*dH);
|
|
|
116 |
|
106 |
Ytemp^[1] := ( - 147.0*Y[0] + 360.0*Y[1] - 450.0*Y[2] + 400.0*Y[3]
|
117 |
Ytemp^[1] := dH60*( - 147.0*Y[0] + 360.0*Y[1] - 450.0*Y[2] + 400.0*Y[3]
|
107 |
- 225.0*Y[4] + 72.0*Y[5] - 10.0*Y[6]) / (60.0*dH);
|
118 |
- 225.0*Y[4] + 72.0*Y[5] - 10.0*Y[6]);
|
108 |
Ytemp^[2] := ( - 10.0*Y[0] - 77.0*Y[1] + 150.0*Y[2] - 100.0*Y[3]
|
119 |
Ytemp^[2] := dH60*( - 10.0*Y[0] - 77.0*Y[1] + 150.0*Y[2] - 100.0*Y[3]
|
109 |
+ 50.0*Y[4] - 15.0*Y[5] + 2.0*Y[6]) / (60.0*dH);
|
120 |
+ 50.0*Y[4] - 15.0*Y[5] + 2.0*Y[6]);
|
110 |
Ytemp^[3] := ( 2.0*Y[0] - 24.0*Y[1] - 35.0*Y[2] + 80.0*Y[3]
|
121 |
Ytemp^[3] := dH60*( + 2.0*Y[0] - 24.0*Y[1] - 35.0*Y[2] + 80.0*Y[3]
|
111 |
- 30.0*Y[4] + 8.0*Y[5] - 1.0*Y[6]) / (60.0*dH);
|
122 |
- 30.0*Y[4] + 8.0*Y[5] - 1.0*Y[6]);
|
112 |
|
123 |
|
113 |
Koeff := (1.0 / (60.0 * dH));
|
124 |
Koeff := (1.0 / (60.0 * dH));
|
114 |
FOR i:=3 TO N-4 DO
|
125 |
FOR i:=3 TO N-4 DO
|
115 |
Ytemp^[i+1] := Koeff*( - 1.0*Y[i-3] + 9.0*Y[i-2]
|
126 |
Ytemp^[i+1] := Koeff*( - 1.0*Y[i-3] + 9.0*Y[i-2]
|
116 |
- 45.0*Y[i-1] + 45.0*Y[i+1]
|
127 |
- 45.0*Y[i-1] + 45.0*Y[i+1]
|
|
... |
|
... |
149 |
8.4887316979573910E+01*Y[N-3] -
|
160 |
8.4887316979573910E+01*Y[N-3] -
|
150 |
6.7909853583659130E+01*Y[N-2] +
|
161 |
6.7909853583659130E+01*Y[N-2] +
|
151 |
2.7729856879994144E+01*Y[N-1]
|
162 |
2.7729856879994144E+01*Y[N-1]
|
152 |
) / (1.1318308930609855E+01*dH);
|
163 |
) / (1.1318308930609855E+01*dH);
|
153 |
FOR i:=1 TO N DO dY1[i-1]:=Ytemp^[i]; END;
|
164 |
FOR i:=1 TO N DO dY1[i-1]:=Ytemp^[i]; END;
|
|
|
165 |
|
154 |
DEALLOCATE(Ytemp,N*TSIZE(LONGREAL));
|
166 |
DEALLOCATE(Ytemp,N*SIZE(LONGREAL));
|
155 |
END Ableit1p7;
|
167 |
END Ableit1p7;
|
156 |
|
168 |
|
157 |
PROCEDURE Ableit2p5(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
169 |
PROCEDURE Ableit2p5(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
158 |
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
|
170 |
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
|
159 |
dH : LONGREAL; (* Schrittweite *)
|
171 |
dH : LONGREAL; (* Schrittweite *)
|
|
... |
|
... |
173 |
IF (n < 5) THEN
|
185 |
IF (n < 5) THEN
|
174 |
Fehler:=TRUE; Fehlerflag:='Zu wenige Datenpunkte (Ableit2p5)';
|
186 |
Fehler:=TRUE; Fehlerflag:='Zu wenige Datenpunkte (Ableit2p5)';
|
175 |
ErrOut(Fehlerflag); RETURN;
|
187 |
ErrOut(Fehlerflag); RETURN;
|
176 |
END;
|
188 |
END;
|
177 |
|
189 |
|
178 |
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
190 |
ALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
179 |
IF (Ytemp = NIL) THEN
|
191 |
IF (Ytemp = NIL) THEN
|
180 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p5) !");
|
192 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p5) !");
|
181 |
END;
|
193 |
END;
|
182 |
|
194 |
|
183 |
recH2 := 1.0 / (dH*dH);
|
195 |
recH2 := 1.0 / (dH*dH);
|
|
... |
|
... |
190 |
END;
|
202 |
END;
|
191 |
Ytemp^[n-1]:= recH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
|
203 |
Ytemp^[n-1]:= recH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
|
192 |
Ytemp^[n ]:= recH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
|
204 |
Ytemp^[n ]:= recH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
|
193 |
|
205 |
|
194 |
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
|
206 |
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
|
|
|
207 |
|
195 |
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
208 |
DEALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
196 |
END Ableit2p5;
|
209 |
END Ableit2p5;
|
197 |
|
210 |
|
198 |
PROCEDURE Ableit2p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
211 |
PROCEDURE Ableit2p7(VAR Y : ARRAY OF LONGREAL; (* ==> Funktionswerte *)
|
199 |
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
|
212 |
VAR dY2 : ARRAY OF LONGREAL; (* <== Zweite Ableitung *)
|
200 |
dH : LONGREAL; (* Schrittweite *)
|
213 |
dH : LONGREAL; (* Schrittweite *)
|
|
... |
|
... |
214 |
END;
|
227 |
END;
|
215 |
IF (n > 8192) THEN
|
228 |
IF (n > 8192) THEN
|
216 |
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit2p7) !");
|
229 |
ErrOut("Zuviele Datenpunkte (n > 8192) (Differ.Ableit2p7) !");
|
217 |
Fehler:=TRUE; RETURN;
|
230 |
Fehler:=TRUE; RETURN;
|
218 |
END;
|
231 |
END;
|
219 |
|
|
|
220 |
ALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
232 |
ALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
221 |
IF (Ytemp = NIL) THEN
|
233 |
IF (Ytemp = NIL) THEN
|
222 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p7) !");
|
234 |
FatalError("Kein Freispeicher vorhanden (Differ.Ableit2p7) !");
|
223 |
END;
|
235 |
END;
|
|
|
236 |
|
224 |
recdH2 := 1.0 / (dH*dH);
|
237 |
recdH2 := 1.0 / (dH*dH);
|
225 |
Ytemp^[1] := recdH2*(2.0*Y[0] - 5.0*Y[1] + 4.0*Y[2] - Y[3]);
|
238 |
Ytemp^[1] := recdH2*(2.0*Y[0] - 5.0*Y[1] + 4.0*Y[2] - Y[3]);
|
226 |
Ytemp^[2] := recdH2*(2.0*Y[1] - 5.0*Y[2] + 4.0*Y[3] - Y[4]);
|
239 |
Ytemp^[2] := recdH2*(2.0*Y[1] - 5.0*Y[2] + 4.0*Y[3] - Y[4]);
|
227 |
rec180dH2 := recdH2/180.0;
|
240 |
rec180dH2 := recdH2/180.0;
|
228 |
FOR i:=2 TO n-5 DO
|
241 |
FOR i:=2 TO n-5 DO
|
|
... |
|
... |
240 |
- Y[n-1]);
|
253 |
- Y[n-1]);
|
241 |
Ytemp^[n-1]:= recdH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
|
254 |
Ytemp^[n-1]:= recdH2*(- Y[n-5] + 4.0*Y[n-4] - 5.0*Y[n-3] + 2.0*Y[n-2]);
|
242 |
Ytemp^[n ]:= recdH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
|
255 |
Ytemp^[n ]:= recdH2*(- Y[n-4] + 4.0*Y[n-3] - 5.0*Y[n-2] + 2.0*Y[n-1]);
|
243 |
|
256 |
|
244 |
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
|
257 |
FOR i:=1 TO n DO dY2[i-1]:=Ytemp^[i]; END;
|
|
|
258 |
|
245 |
DEALLOCATE(Ytemp,n*TSIZE(LONGREAL));
|
259 |
DEALLOCATE(Ytemp,n*SIZE(LONGREAL));
|
246 |
END Ableit2p7;
|
260 |
END Ableit2p7;
|
247 |
|
261 |
|
248 |
PROCEDURE AbleitN(VAR X,Y : ARRAY OF LONGREAL;
|
262 |
PROCEDURE AbleitN(VAR X,Y : ARRAY OF LONGREAL;
|
249 |
VAR dY : ARRAY OF LONGREAL;
|
263 |
VAR dY : ARRAY OF LONGREAL;
|
250 |
ndata : CARDINAL;
|
264 |
ndata : CARDINAL;
|
|
... |
|
... |
266 |
IF (Grad < NAbl+1) THEN Grad:=NAbl+1; END;
|
280 |
IF (Grad < NAbl+1) THEN Grad:=NAbl+1; END;
|
267 |
dim := Grad + 2;
|
281 |
dim := Grad + 2;
|
268 |
IF (ndata < dim) THEN
|
282 |
IF (ndata < dim) THEN
|
269 |
FatalError("Zuwenig Datenpunkte (Ableit1) !");
|
283 |
FatalError("Zuwenig Datenpunkte (Ableit1) !");
|
270 |
END;
|
284 |
END;
|
271 |
ALLOCATE(dYloc ,ndata *TSIZE(LONGREAL));
|
285 |
ALLOCATE(dYloc ,ndata *SIZE(LONGREAL));
|
272 |
ALLOCATE(PKoeff,(Grad+1)*TSIZE(LONGREAL));
|
286 |
ALLOCATE(PKoeff,(Grad+1)*SIZE(LONGREAL));
|
273 |
IF (PKoeff = NIL) OR (dYloc = NIL) THEN
|
287 |
IF (PKoeff = NIL) OR (dYloc = NIL) THEN
|
274 |
FatalError("Kein Freispeicher vorhanden (Ableit1) !");
|
288 |
FatalError("Kein Freispeicher vorhanden (Ableit1) !");
|
275 |
END;
|
289 |
END;
|
276 |
|
290 |
|
277 |
grad := Grad; (* wg. Compilerwarnung *)
|
291 |
grad := Grad; (* wg. Compilerwarnung *)
|
|
... |
|
... |
318 |
|
332 |
|
319 |
FOR i:=0 TO ndata-1 DO (* Kopiere L"osung *)
|
333 |
FOR i:=0 TO ndata-1 DO (* Kopiere L"osung *)
|
320 |
dY[i]:=dYloc^[i];
|
334 |
dY[i]:=dYloc^[i];
|
321 |
END;
|
335 |
END;
|
322 |
|
336 |
|
323 |
DEALLOCATE(dYloc ,ndata *TSIZE(LONGREAL));
|
337 |
DEALLOCATE(dYloc ,ndata *SIZE(LONGREAL));
|
324 |
DEALLOCATE(PKoeff,(Grad+1)*TSIZE(LONGREAL));
|
338 |
DEALLOCATE(PKoeff,(Grad+1)*SIZE(LONGREAL));
|
325 |
END AbleitN;
|
339 |
END AbleitN;
|
326 |
|
340 |
|
327 |
PROCEDURE Ableit1und2( X : LONGREAL;
|
341 |
PROCEDURE Ableit1und2( X : LONGREAL;
|
328 |
H : LONGREAL;
|
342 |
H : LONGREAL;
|
329 |
Funk : Funktion1;
|
343 |
Funk : Funktion1;
|