Integral.mod.m2pp
539 lines (483 with data), 18.4 kB
IMPLEMENTATION MODULE Integral;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei Integral.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < Integral.mod.m2pp > Integral.mod *)
(* *)
(* mit Parameter = {XDS|GM2} erzeugt werden koennen. *)
(* *)
(* XDS : Der Export von Objekten aus lokalen Modulen wird mit *)
(* "IMPORT" angezeigt - das scheint nicht ganz "regelkonform" *)
(* GM2 : Normales "EXPORT" wird genutzt (ELSE-Fall). *)
(* *)
(* ansonsten gibt es keine Aenderungen am Quellcode *)
(* *)
(* Fuer die Routine InitGausInt kann Debug-Code mit der Option *)
(* -D __DEBUG__ eingefuegt werden *)
(*------------------------------------------------------------------------*)
(* Module providing procedures for numerical integration *)
(*------------------------------------------------------------------------*)
(* Letzt Bearbeitung: *)
(* *)
(* 17.07.93, MRi: Durchsicht *)
(* 01.03.18, MRi: Kosmetische Korrekturen, TriQR durch TriQL ersetzt *)
(* 02.03.18, MRi: Korrektur in NumInt, Simson & SimCum eingefuegt *)
(*------------------------------------------------------------------------*)
(* Testprogram in TstIntegral *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: Integral.mod.m2pp,v 1.5 2018/03/02 10:48:20 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps,Funktion1;
FROM Deklera IMPORT VEKTOR,CARDVEKTOR,MaxDim;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
FROM MatLib IMPORT QuickSort;
FROM EigenLib1 IMPORT ListenZeiger,TriQL,BerechneEV;
<* IF (__DEBUG__) THEN *>
IMPORT TIO;
<* END *>
MODULE GAUSSINT;
IMPORT VEKTOR,CARDVEKTOR,MaxDim;
IMPORT QuickSort;
IMPORT ListenZeiger,TriQL,BerechneEV;
IMPORT sqrt;
IMPORT MachEps,Funktion1;
<* IF (__DEBUG__) THEN *>
IMPORT TIO;
<* END *>
<* IF (__XDS__) THEN *>
IMPORT InitGaussInt,GaussInt; (* EXPORT *)
<* ELSE *>
EXPORT InitGaussInt,GaussInt;
<* END *>
(*------------------ Lokale Objekte ------------------------------*)
CONST genau = 10.0*MachEps;
VAR A,X : VEKTOR;
Anzahl : CARDINAL; (* Anzahl der Stuetzstellen *)
(*----------------------------------------------------------------*)
PROCEDURE InitGaussInt(NSt : CARDINAL);
VAR HD,ND,EVi : VEKTOR;
I : LONGREAL;
i : CARDINAL;
ListPtr : ListenZeiger;
II : CARDVEKTOR;
BEGIN
IF (NSt <= 1) THEN X[1]:=0.0; A[1]:=2.0; Anzahl:=1; RETURN; END;
IF (NSt > MaxDim-2) THEN NSt:=MaxDim-2; END;
Anzahl:=NSt;
FOR i:=1 TO NSt-1 DO
I:=VAL(LONGREAL,i);
ND[i]:=I / sqrt(4.0*I*I - 1.0);
HD[i]:=0.0;
END;
ND[NSt]:=0.0;
HD[NSt]:=0.0;
(* Die X[i] sind die Eigenwerte der Trilinearmatrix *)
TriQL(HD,ND,X,NSt,ListPtr,genau);
QuickSort(X,II,NSt,1);
FOR i:=1 TO NSt DO
BerechneEV(HD,ND,X[i],EVi,II[i],genau,30,ListPtr,NSt);
A[i]:=2.0*EVi[1]*EVi[1];
<* IF (__DEBUG__) THEN *>
TIO.WrLngReal(X[i],24,-16);
TIO.WrLngReal(A[i],24,-16);
TIO.WrLn;
<* END *>
END;
END InitGaussInt;
PROCEDURE GaussInt(Untergr,Obergr : LONGREAL;
N : CARDINAL;
Funk : Funktion1) : LONGREAL;
VAR Sum : LONGREAL;
p,h : LONGREAL;
i,k : CARDINAL;
BEGIN
h:=ABS(Obergr-Untergr) / VAL(LONGREAL,2*N);
Sum:=0.0; p:=Untergr + h;
FOR i:=1 TO N DO
FOR k:=1 TO Anzahl DO Sum:=Sum + A[k]*h*Funk(p + X[k]*h); END;
p:=p + 2.0*h;
END;
RETURN Sum;
END GaussInt;
BEGIN (* Zur Sicherheit vorbelegt *)
Anzahl:=2;
X[1]:=0.577350269189626E+00; X[2]:=-X[1];
A[1]:=1.000000000000000E+00; A[2]:=A[1];
END GAUSSINT;
PROCEDURE GaussIntKonst(Untergr,Obergr : LONGREAL;
N : CARDINAL;
NSt : CARDINAL; (* Ordnung der Integration *)
Funk : Funktion1) : LONGREAL;
VAR Sum : LONGREAL;
p,h : LONGREAL;
A,X : ARRAY [1..24] OF LONGREAL;
i,k : CARDINAL;
BEGIN
IF NOT (NSt IN {4,6,8,10,12,16,24}) THEN
ErrOut("Anzahl Stützpunkte nicht implementiert, nehme 4");
NSt:=4;
END;
CASE NSt OF
2: (* GAUSS-02 DAT *)
X[1]:=0.577350269189626E+00; X[2]:=-X[1];
A[1]:=1.000000000000000E+00; A[2]:=A[1];
| 4: (* GAUSS-04 DAT *)
X[1]:=0.861136311594053E+00; X[4]:=-X[1];
X[2]:=0.339981043584856E+00; X[3]:=-X[2];
A[1]:=0.347854845137454E+00; A[4]:=A[1];
A[2]:=0.652145154862546E+00; A[3]:=A[2];
| 6: (* GAUSS-06 DAT *)
A[1]:=0.171324492379170E+00; A[6]:=A[1];
A[2]:=0.360761573048139E+00; A[5]:=A[2];
A[3]:=0.467913934572691E+00; A[4]:=A[3];
X[1]:=0.932469514203152E+00; X[6]:=-X[1];
X[2]:=0.661209386466265E+00; X[5]:=-X[2];
X[3]:=0.238619186083197E+00; X[4]:=-X[3];
| 8: (* GAUSS-08 DAT *)
X[1]:=0.960289856497536E+00; X[8]:=-X[1];
X[2]:=0.796666477413627E+00; X[7]:=-X[2];
X[3]:=0.525532409916329E+00; X[6]:=-X[3];
X[4]:=0.183434642495650E+00; X[5]:=-X[4];
A[1]:=0.101228536290376E+00; A[8]:=A[1];
A[2]:=0.222381034453374E+00; A[7]:=A[2];
A[3]:=0.313706645877887E+00; A[6]:=A[3];
A[4]:=0.362683783378362E+00; A[5]:=A[4];
|10: (* GAUSS-10 DAT *)
X[1]:=0.973906528517172E+00; X[10]:=-X[1];
X[2]:=0.865063366688985E+00; X[ 9]:=-X[2];
X[3]:=0.679409568299024E+00; X[ 8]:=-X[3];
X[4]:=0.433395394129247E+00; X[ 7]:=-X[4];
X[5]:=0.148874338981631E+00; X[ 6]:=-X[5];
A[1]:=0.666713443086880E-01; A[10]:=A[1];
A[2]:=0.149451349150581E+00; A[ 9]:=A[2];
A[3]:=0.219086362515982E+00; A[ 8]:=A[3];
A[4]:=0.269266719309996E+00; A[ 7]:=A[4];
A[5]:=0.295524224714753E+00; A[ 6]:=A[5];
|12: (* GAUSS-12 DAT *)
X[1]:=0.981560634246719E+00; X[12]:=-X[1];
X[2]:=0.904117256370475E+00; X[11]:=-X[2];
X[3]:=0.769902674194305E+00; X[10]:=-X[3];
X[4]:=0.587317954286617E+00; X[ 9]:=-X[4];
X[5]:=0.367831498998180E+00; X[ 8]:=-X[5];
X[6]:=0.125233408511469E+00; X[ 7]:=-X[6];
A[1]:=0.471753363865120E-01; A[12]:=A[1];
A[2]:=0.106939325995318E+00; A[11]:=A[2];
A[3]:=0.160078328543346E+00; A[10]:=A[3];
A[4]:=0.203167426723066E+00; A[ 9]:=A[4];
A[5]:=0.233492536538355E+00; A[ 8]:=A[5];
A[6]:=0.249147045813403E+00; A[ 7]:=A[6];
| 16: (* GAUSS-16 DAT *)
X[1]:=0.989400934991649939E+00; X[16]:=-X[1];
X[2]:=0.944575023073232573E+00; X[15]:=-X[2];
X[3]:=0.865631202387831741E+00; X[14]:=-X[3];
X[4]:=0.755404408355003026E+00; X[13]:=-X[4];
X[5]:=0.617876244402643743E+00; X[12]:=-X[5];
X[6]:=0.458016777657227384E+00; X[11]:=-X[6];
X[7]:=0.281603550779258915E+00; X[10]:=-X[7];
X[8]:=0.950125098376374405E-01; X[ 9]:=-X[8];
A[1]:=0.271524594117540947E-01; A[16]:=A[1];
A[2]:=0.622535239386478928E-01; A[15]:=A[2];
A[3]:=0.951585116824927857E-01; A[14]:=A[3];
A[4]:=0.124628971255533877E+00; A[13]:=A[4];
A[5]:=0.149595988816576736E+00; A[12]:=A[5];
A[6]:=0.169156519395002536E+00; A[11]:=A[6];
A[7]:=0.182603415044923584E+00; A[10]:=A[7];
A[8]:=0.189450610455068502E+00; A[ 9]:=A[8];
|24: (* GAUSS-24 DAT *)
A[ 1]:=0.123412297999871993E-01; A[24]:=A[1];
A[ 2]:=0.285313886289336634E-01; A[23]:=A[2];
A[ 3]:=0.442774388174198060E-01; A[22]:=A[3];
A[ 4]:=0.592985849154367807E-01; A[21]:=A[4];
A[ 5]:=0.733464814110802998E-01; A[20]:=A[5];
A[ 6]:=0.861901615319532743E-01; A[19]:=A[6];
A[ 7]:=0.976186521041138844E-01; A[18]:=A[7];
A[ 8]:=0.107444270115965634E+00; A[17]:=A[8];
A[ 9]:=0.115505668053725599E+00; A[16]:=A[9];
A[10]:=0.121670472927803391E+00; A[15]:=A[10];
A[11]:=0.125837456346828303E+00; A[14]:=A[11];
A[12]:=0.127938195346752159E+00; A[13]:=A[12];
X[ 1]:=0.995187219997021366E+00; X[24]:=-X[1];
X[ 2]:=0.974728555971309502E+00; X[23]:=-X[2];
X[ 3]:=0.938274552002732756E+00; X[22]:=-X[3];
X[ 4]:=0.886415527004401030E+00; X[21]:=-X[4];
X[ 5]:=0.820001985973902919E+00; X[20]:=-X[5];
X[ 6]:=0.740124191578554358E+00; X[19]:=-X[6];
X[ 7]:=0.648093651936975573E+00; X[18]:=-X[7];
X[ 8]:=0.545421471388839535E+00; X[17]:=-X[8];
X[ 9]:=0.433793507626045141E+00; X[16]:=-X[9];
X[10]:=0.315042679696163369E+00; X[15]:=-X[10];
X[11]:=0.191118867473616311E+00; X[14]:=-X[11];
X[12]:=0.640568928626056300E-01; X[13]:=-X[12];
ELSE
END;
h:=ABS(Obergr-Untergr) / VAL(LONGREAL,2*N);
Sum:=0.0; p:=Untergr + h;
FOR i:=1 TO N DO
FOR k:=1 TO NSt DO Sum:=Sum + A[k]*h*Funk(p + X[k]*h); END;
p:=p + 2.0*h;
END;
RETURN Sum;
END GaussIntKonst;
PROCEDURE SimpInt(Untergr,Obergr : LONGREAL;
genau : LONGREAL;
MaxIter : CARDINAL;
VAR Iter : CARDINAL;
Funk : Funktion1) : LONGREAL;
VAR i,n : CARDINAL;
Sum,ss,S1,S2,S4, H,H2,x : LONGREAL;
BEGIN
n:=2; Fehler:=FALSE;
H:=(Obergr-Untergr)*0.5;
S1:=H*(Funk(Untergr) + Funk(Obergr));
S2:=0.0;
S4:=4.0*H*Funk(Untergr+H);
Sum:=S1*S2*S4;
Iter:=0;
REPEAT
INC(Iter);
ss:=Sum;
n:=2*n;
S2:=0.5*S2 + 0.25*S4;
S1:=0.5*S1; H:=0.5*H;
i:=1; S4:=0.0;
H2:=2.0*H; x:=Untergr+H;
REPEAT
S4:=S4 + Funk(x);
x:=x+H2;
INC(i,2);
UNTIL (i > n);
S4:=4.0*H*S4;
Sum:=S1 + S2 + S4;
UNTIL (ABS(Sum-ss) < genau*ABS(Sum)) OR (Iter >= MaxIter);
IF (Iter >= MaxIter) THEN
Fehler:=TRUE;
Fehlerflag:='MaxIter "uberschritten (SimpInt)';
ErrOut(Fehlerflag);
END;
RETURN Sum / 3.0;
END SimpInt;
PROCEDURE Romberg(UnterGr : LONGREAL;
OberGr : LONGREAL;
genau : LONGREAL;
MaxIter : CARDINAL;
VAR Iter : CARDINAL;
Funk : Funktion1) : LONGREAL;
CONST eps = 1.0E-12;
Max = 100;
TYPE Tableau = ARRAY [1..Max] OF LONGREAL;
VAR Xneu,Xalt : Tableau; (* Integrationstableaus *)
i,j,ProdIter : CARDINAL;
dH,Sum,Prod : LONGREAL; (* dH : Abstand d. Int.-punkte *)
BEGIN
IF (genau <= eps) THEN genau:=eps; END;
IF (MaxIter = 0) THEN MaxIter:=100; END;
IF (MaxIter > Max) THEN MaxIter:=Max; END;
dH:=OberGr - UnterGr;
Xalt[1]:=0.5*dH*(Funk(UnterGr) + Funk(OberGr));
Iter:=1; ProdIter:=1;
REPEAT
INC(Iter);
Sum:=0.0;
FOR i:=1 TO ProdIter DO
Sum:=Sum + Funk(UnterGr + (VAL(LONGREAL,i) - 0.5)*dH);
END;
Xneu[1]:=0.5*(Xalt[1] + dH*Sum);
ProdIter:=2*ProdIter;
Prod:=1.0;
FOR j:=2 TO Iter DO
Prod:=4.0*Prod;
Xneu[j]:=(Prod*Xneu[j-1] - Xalt[j-1]) / (Prod - 1.0);
END;
dH:=0.5*dH;
Xalt:=Xneu;
UNTIL (ABS(Xneu[Iter - 1] - Xneu[Iter]) <= ABS(genau * Xneu[Iter]))
OR (Iter >= MaxIter);
IF Iter >= MaxIter THEN
Fehler:=TRUE; Fehlerflag:='MaxIter "uberschritten (Romberg).';
ErrOut(Fehlerflag);
END;
RETURN Xneu[Iter];
END Romberg;
PROCEDURE Quadrature(Untergr,Obergr : LONGREAL;
eps : LONGREAL; (* Geforderte Genauigkeit *)
VAR err : LONGREAL; (* Fehlerabsch"atzung *)
Funkt : Funktion1) : LONGREAL;
CONST eta = 1.0E-15;
VAR H,c,d1,den,ddt,e,gr,hm,
nt,sm,t,t1,t2,t2a,ta,tab,tb,v,w : LONGREAL;
int : LONGREAL;
dt : ARRAY [0..6] OF LONGREAL;
d : ARRAY [1..6] OF LONGREAL;
i,m,mr,n,nn : CARDINAL;
bo,bu,odd : BOOLEAN;
BEGIN
Fehler:=FALSE;
IF (eps < eta) THEN eps:=eta; END;
FOR i:=0 TO 6 DO dt[i]:=0.0; END;
FOR i:=1 TO 6 DO d[i]:=0.0; END;
n:=2; nn:=3;
H:=Obergr-Untergr;
t1:=0.0; gr:=0.0; sm:=0.0;
t2:=0.5*(Funkt(Untergr) + Funkt(Obergr)); t2a:=t2; tb:=ABS(t2);
c:=t2*H; dt[0]:=c;
odd:=TRUE; bu:=FALSE;
m:=0;
LOOP
INC(m);
bo:= (m >= 7);
hm:=H/VAL(LONGREAL,n);
IF odd THEN
FOR i:=1 TO n BY 2 DO
w:=Funkt(Untergr+VAL(LONGREAL,i)*hm); (**)
t2:=t2+w;
tb:=tb+ABS(w);
END;
nt:=t2;
tab:=tb*ABS(hm);
d[1]:= 16.0/9.0;
d[3]:= 64.0/9.0;
d[5]:=256.0/9.0;
ELSE
FOR i:=1 TO n BY 6 DO
w:=VAL(LONGREAL,i)*hm;
t1:=Funkt(Untergr+w) + Funkt(Obergr-w) + t1;
END;
nt:=t1+t2a;
t2a:=t2;
d[1]:=9.0/4.0;
d[3]:=9.0;
d[5]:=36.0;
END; (* IF odd *)
ddt:=dt[0];
nt:=nt*hm; t:=nt; dt[0]:=nt;
IF bo THEN
mr:=6;
d[6]:=64.0;
w:=144.0;
ELSE
mr:=m;
w:=VAL(LONGREAL,n*n); d[m]:=w;
END; (* IF bo *)
FOR i:=1 TO mr DO
d1:=d[i]*ddt;
den:=d1-nt;
e:=nt-ddt;
IF (den <> 0.0) THEN
e:=e/den;
v:=nt*e;
nt:=d1*e;
t:=t+v;
ELSE
v:=0.0;
nt:=0.0;
END; (* IF *)
ddt:=dt[i];
dt[i]:=v;
END; (* FOR i *)
ta:=c;
int:=t; c:=t;
IF NOT bo THEN t:=t-v; END;
v:=t-ta;
t:=t+v;
err:=ABS(v);
IF (ta < t) THEN
d1:=ta;
ta:=t;
t:=d1;
END;
bo:=bo AND ((ta < gr) OR (t > sm));
IF bu OR bo OR (err < eps*tab*w) THEN EXIT END;
gr:=ta;
sm:=t;
odd:=NOT odd;
i:=n;
n:=nn;
nn:=i*i;
bu:=bo;
d[2]:=4.0;
d[4]:=16.0;
IF (m = 15) THEN
Fehler:=TRUE;
Fehlerflag:='keine Konvergenz (Quadrature)';
ErrOut(Fehlerflag);
EXIT;
END;
END; (* LOOP m *)
v:=eta*tab;
IF (err < v) THEN err:=v; END;
IF (ABS(err) > ABS(eps)) AND NOT Fehler THEN
Fehler:=TRUE;
Fehlerflag:='|err| > |eps| (Quadrature)';
ErrOut(Fehlerflag);
END;
RETURN int;
END Quadrature;
PROCEDURE NumInt(VAR X,Y : ARRAY OF LONGREAL;
N : CARDINAL) : LONGREAL;
VAR i : CARDINAL;
dX,Sum : LONGREAL;
BEGIN
Sum:=0.0;
FOR i:=1 TO N-1 DO
dX := 0.5*(X[i] -X[i-1]);
Sum:=Sum + dX*(Y[i] + Y[i-1]);
END;
RETURN Sum;
END NumInt;
PROCEDURE Simson(VAR Y : ARRAY OF LONGREAL;
H : LONGREAL;
N : CARDINAL) : LONGREAL;
VAR i,j : CARDINAL;
sum : LONGREAL;
BEGIN
sum:=0.0;
j:=0;
FOR i:=2 TO N-1 BY 2 DO
sum:=sum + Y[j] + 4.0*Y[i-1] + Y[i];
j:=i;
END;
IF NOT ODD(N) THEN
sum:=sum - 0.25*Y[N-3] + 2.0*Y[j] + 1.25*Y[N-1];
END;
RETURN sum*H / 3.0;
END Simson;
PROCEDURE SimCum(VAR Y : ARRAY OF LONGREAL;
H : LONGREAL;
N : CARDINAL) : LONGREAL;
VAR i,j,k : CARDINAL;
fac,S : LONGREAL;
term : ARRAY [1..5] OF LONGREAL;
Area : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
BEGIN
ALLOCATE(Area,N*TSIZE(LONGREAL));
Area^[0] := 0.0;
fac := H / 3.0;
term[1] := 1.25*Y[0];
term[2] := - 0.25*Y[0];
k:=0;
FOR i:=2 TO N-1 BY 2 DO
j := i - 1;
term[3] := 2.00*Y[j];
term[4] := 1.25*Y[i];
term[5] := - 0.25*Y[i];
Area^[j] := Area^[k] + fac*(term[1] + term[3] + term[5]);
Area^[i] := Area^[j] + fac*(term[2] + term[3] + term[4]);
term[1] := term[4];
term[2] := term[5];
k := i;
END; (* FOR *)
IF NOT ODD(N) THEN
Area^[N-1]:=Area^[k] + fac*(-0.25*Y[j] + 2.0*Y[k] + 1.25*Y[N-1]);
END;
S := Area^[N-1];
DEALLOCATE(Area,N*TSIZE(LONGREAL));
RETURN S;
END SimCum;
END Integral.