MODULE CalcGauss;
IMPORT NbrInt, NbrRe, NbrCplx, Data, DataLists, CalcFn;
CONST
Coarse* = 99; Medium* = 100; Fine* = 101;
OKay* = 0; MaxSubDivReached* = 1; RoundoffError* = 2; RoughIntegrand* = 3;
VAR
MaxIntervals-: NbrInt.Integer;
node8, wgtGauss4, wgtKronrod8, node16, wgtGauss8, wgtKronrod16, node31, wgtGauss15,
wgtKronrod31: POINTER TO ARRAY OF NbrRe.Real;
TYPE
Interval = OBJECT (Data.Datum)
VAR a, b, error: NbrRe.Real;
PROCEDURE & Initialize*;
BEGIN
Initialize^;
a := 0; b := 0; error := 0
END Initialize;
PROCEDURE Copy( VAR copy: Data.Datum );
VAR new, obj: Interval;
BEGIN
IF (copy = NIL ) OR ~(copy IS Interval) THEN NEW( new ); copy := new END;
Copy^( copy );
obj := copy( Interval );
obj.a := a; obj.b := b; obj.error := error;
copy := obj
END Copy;
END Interval;
ReInterval = OBJECT (Interval)
VAR soln: NbrRe.Real;
PROCEDURE & Initialize*;
BEGIN
Initialize^;
soln := 0
END Initialize;
PROCEDURE Copy( VAR copy: Data.Datum );
VAR new, obj: ReInterval;
BEGIN
IF (copy = NIL ) OR ~(copy IS ReInterval) THEN NEW( new ); copy := new END;
Copy^( copy );
obj := copy( ReInterval );
obj.soln := soln;
copy := obj
END Copy;
END ReInterval;
CplxInterval = OBJECT (Interval)
VAR soln: NbrCplx.Complex;
PROCEDURE & Initialize*;
BEGIN
Initialize^;
soln := 0
END Initialize;
PROCEDURE Copy( VAR copy: Data.Datum );
VAR new, obj: CplxInterval;
BEGIN
IF (copy = NIL ) OR ~(copy IS CplxInterval) THEN NEW( new ); copy := new END;
Copy^( copy );
obj := copy( CplxInterval );
obj.soln := soln;
copy := obj
END Copy;
END CplxInterval;
PROCEDURE GetKey( a, b, atX: NbrRe.Real ): Data.Key;
VAR k: Data.Key;
BEGIN
k := NbrRe.Int( MaxIntervals * (atX - a)/(b - a) ); RETURN k
END GetKey;
PROCEDURE ReGaussKronrod( f: CalcFn.ReArg; fromX, toX: NbrRe.Real; integrator: NbrInt.Integer;
VAR result, absError, absResult: NbrRe.Real );
VAR i, len: NbrInt.Integer;
abscissa, center, fCenter, fSum, halfLength, resultGauss, resultKronrod, tolerance: NbrRe.Real;
fAbove, fBelow, node, wgtGauss, wgtKronrod: POINTER TO ARRAY OF NbrRe.Real;
BEGIN
IF integrator = Coarse THEN
len := 8; NEW( node, 8 ); NEW( wgtGauss, 4 ); NEW( wgtKronrod, 8 );
FOR i := 0 TO 7 DO node[i] := node8[i]; wgtKronrod[i] := wgtKronrod8[i] END;
FOR i := 0 TO 3 DO wgtGauss[i] := wgtGauss4[i] END
ELSIF integrator = Fine THEN
len := 31; NEW( node, 31 ); NEW( wgtGauss, 15 ); NEW( wgtKronrod, 31 );
FOR i := 0 TO 30 DO node[i] := node31[i]; wgtKronrod[i] := wgtKronrod31[i] END;
FOR i := 0 TO 14 DO wgtGauss[i] := wgtGauss15[i] END
ELSE
len := 16; NEW( node, 16 ); NEW( wgtGauss, 8 ); NEW( wgtKronrod, 16 );
FOR i := 0 TO 15 DO node[i] := node16[i]; wgtKronrod[i] := wgtKronrod16[i] END;
FOR i := 0 TO 7 DO wgtGauss[i] := wgtGauss8[i] END
END;
NEW( fAbove, len - 1 ); NEW( fBelow, len - 1 ); center := (fromX + toX) / 2; halfLength := (toX - fromX) / 2;
fCenter := f( center ); resultKronrod := wgtKronrod[len - 1] * fCenter; absResult := ABS( resultKronrod );
IF NbrInt.Odd( len ) THEN resultGauss := 0 ELSE resultGauss := wgtGauss[(len DIV 2) - 1] * fCenter END;
FOR i := 0 TO len - 2 DO
abscissa := halfLength * node[i]; fAbove[i] := f( center + abscissa ); fBelow[i] := f( center - abscissa );
fSum := fAbove[i] + fBelow[i]; resultKronrod := resultKronrod + wgtKronrod[i] * fSum;
absResult := absResult + wgtKronrod[i] * (ABS( fAbove[i] ) + ABS( fBelow[i] ));
IF NbrInt.Odd( i ) THEN resultGauss := resultGauss + wgtGauss[i DIV 2] * fSum END
END;
result := halfLength * resultKronrod; absResult := halfLength * absResult;
absError := ABS( halfLength * (resultKronrod - resultGauss) ); tolerance := 50 * NbrRe.Epsilon;
IF absResult > 1/NbrRe.MaxNbr*tolerance THEN absError := NbrRe.Max( absError, tolerance * absResult ) END
END ReGaussKronrod;
PROCEDURE Solve*( f: CalcFn.ReArg; a, b: NbrRe.Real; integrator: NbrInt.Integer; VAR error: NbrRe.Real;
VAR res: NbrInt.Integer ): NbrRe.Real;
VAR ignor, successful: BOOLEAN; key, maxKey: Data.Key; sign: NbrInt.Integer;
aa, bb, maxError, maxTol, midPoint, sumError, sumResult, tolerance: NbrRe.Real; datum: Data.Datum;
interval, intervalL, intervalR: ReInterval; history: DataLists.List;
PROCEDURE Create( withKey: Data.Key ): ReInterval;
VAR int: ReInterval;
BEGIN
NEW( int ); int.SetKey( withKey ); RETURN int
END Create;
PROCEDURE Update( VAR int: ReInterval );
VAR absApprox, absError, approx: NbrRe.Real;
BEGIN
ReGaussKronrod( f, int.a, int.b, integrator, approx, absError, absApprox ); int.soln := approx;
int.error := absError * (int.b - int.a) / (bb - aa);
IF (absError < 100 * NbrRe.Epsilon * absApprox) & (absError > tolerance) THEN
res := RoundoffError
END;
IF NbrRe.Abs( int.b ) < ((1 + 100*NbrRe.Epsilon) * (NbrRe.Abs( int.a ) + 1000/NbrRe.MaxNbr)) THEN
res := RoughIntegrand
END
END Update;
BEGIN
IF f # NIL THEN
res := OKay;
IF a = b THEN error := 0; sign := 1; sumResult := 0
ELSE
maxTol := 0.1; tolerance := NbrRe.Max( 100 * NbrRe.Epsilon, NbrRe.Min( maxTol, error ) );
IF a < b THEN sign := 1; aa := a; bb := b ELSE sign := -1; aa := b; bb := a END;
key := GetKey( aa, bb, bb ); interval := Create( key ); interval.a := aa; interval.b := bb;
Update( interval ); sumResult := interval.soln; sumError := interval.error;
IF NbrRe.Abs( sumResult ) < 1 THEN error := sumError
ELSE error := NbrRe.Abs( sumError / sumResult )
END;
IF error > tolerance THEN
NEW( history ); history.Insert( interval, ignor );
LOOP
history.rider.Home; datum := history.rider.Inspect(); interval := datum( ReInterval );
maxError := interval.error; interval.GetKey( maxKey );
WHILE ~history.rider.eol DO
history.rider.Next; datum := history.rider.Inspect(); interval := datum( ReInterval );
IF interval.error > maxError THEN
maxError := interval.error; interval.GetKey( maxKey )
END
END;
datum := history.rider.Retrieve( maxKey ); interval := datum( ReInterval );
history.Delete( maxKey, ignor ); midPoint := (interval.a + interval.b) / 2;
key := GetKey( interval.a, interval.b, midPoint ); intervalL := Create( key );
intervalL.a := interval.a; intervalL.b := midPoint; Update( intervalL );
history.Insert( intervalL, successful );
IF successful & (key # maxKey) THEN
intervalR := Create( maxKey ); intervalR.a := midPoint; intervalR.b := interval.b;
Update( intervalR ); history.Insert( intervalR, ignor )
ELSE
IF successful THEN history.Delete( key, ignor ) END;
history.Insert( interval, ignor ); res := MaxSubDivReached
END;
IF res # OKay THEN EXIT END;
sumResult := sumResult + intervalL.soln + intervalR.soln - interval.soln;
sumError := sumError + intervalL.error + intervalR.error - interval.error;
IF NbrRe.Abs( sumResult ) < 1 THEN error := sumError
ELSE error := NbrRe.Abs( sumError / sumResult )
END;
IF error < tolerance THEN EXIT END
END
END
END
ELSE sign := 1; sumResult := 0
END;
RETURN sign * sumResult
END Solve;
PROCEDURE CplxGaussKronrod( f: CalcFn.MixedArg; fromX, toX: NbrRe.Real; z: NbrCplx.Complex;
integrator: NbrInt.Integer; VAR result: NbrCplx.Complex;
VAR absError, absResult: NbrRe.Real );
VAR i, len: NbrInt.Integer; abscissa, center, halfLength, tolerance: NbrRe.Real;
fCenter, fSum, resultGauss, resultKronrod: NbrCplx.Complex;
node, wgtGauss, wgtKronrod: POINTER TO ARRAY OF NbrRe.Real;
fAbove, fBelow: POINTER TO ARRAY OF NbrCplx.Complex;
BEGIN
IF integrator = Coarse THEN
len := 8; NEW( node, 8 ); NEW( wgtGauss, 4 ); NEW( wgtKronrod, 8 );
FOR i := 0 TO 7 DO node[i] := node8[i]; wgtKronrod[i] := wgtKronrod8[i] END;
FOR i := 0 TO 3 DO wgtGauss[i] := wgtGauss4[i] END
ELSIF integrator = Fine THEN
len := 31; NEW( node, 31 ); NEW( wgtGauss, 15 ); NEW( wgtKronrod, 31 );
FOR i := 0 TO 30 DO node[i] := node31[i]; wgtKronrod[i] := wgtKronrod31[i] END;
FOR i := 0 TO 14 DO wgtGauss[i] := wgtGauss15[i] END
ELSE
len := 16; NEW( node, 16 ); NEW( wgtGauss, 8 ); NEW( wgtKronrod, 16 );
FOR i := 0 TO 15 DO node[i] := node16[i]; wgtKronrod[i] := wgtKronrod16[i] END;
FOR i := 0 TO 7 DO wgtGauss[i] := wgtGauss8[i] END
END;
NEW( fAbove, len - 1 ); NEW( fBelow, len - 1 ); center := (fromX + toX) / 2; halfLength := (toX - fromX) / 2;
fCenter := f( center, z ); resultKronrod := wgtKronrod[len - 1] * fCenter;
absResult := NbrCplx.Abs( resultKronrod );
IF NbrInt.Odd( len ) THEN resultGauss := 0 ELSE resultGauss := wgtGauss[(len DIV 2) - 1] * fCenter END;
FOR i := 0 TO len - 2 DO
abscissa := halfLength * node[i]; fAbove[i] := f( center + abscissa, z ); fBelow[i] := f( center - abscissa, z );
fSum := fAbove[i] + fBelow[i]; resultKronrod := resultKronrod + wgtKronrod[i] * fSum;
absResult := absResult + wgtKronrod[i] * (NbrCplx.Abs( fAbove[i] ) + NbrCplx.Abs( fBelow[i] ));
IF NbrInt.Odd( i ) THEN resultGauss := resultGauss + wgtGauss[i DIV 2] * fSum END
END;
result := halfLength * resultKronrod; absResult := halfLength * absResult;
absError := NbrCplx.Abs( halfLength * (resultKronrod - resultGauss) ); tolerance := 50 * NbrRe.Epsilon;
IF absResult > 1/NbrRe.MaxNbr*tolerance THEN absError := NbrRe.Max( absError, tolerance * absResult ) END
END CplxGaussKronrod;
PROCEDURE SolveCplx*( f: CalcFn.MixedArg; a, b: NbrRe.Real; z: NbrCplx.Complex; integrator: NbrInt.Integer;
VAR error: NbrRe.Real; VAR res: NbrInt.Integer ): NbrCplx.Complex;
VAR ignor, successful: BOOLEAN; key, maxKey: Data.Key; sign: NbrInt.Integer;
aa, bb, maxError, maxTol, midPoint, sumError, tolerance: NbrRe.Real; sumResult: NbrCplx.Complex;
datum: Data.Datum; interval, intervalL, intervalR: CplxInterval; history: DataLists.List;
PROCEDURE Create( withKey: Data.Key ): CplxInterval;
VAR int: CplxInterval;
BEGIN
NEW( int ); int.SetKey( withKey ); RETURN int
END Create;
PROCEDURE Update( VAR int: CplxInterval );
VAR absApprox, absError: NbrRe.Real; approx: NbrCplx.Complex;
BEGIN
CplxGaussKronrod( f, int.a, int.b, z, integrator, approx, absError, absApprox ); int.soln := approx;
int.error := absError * (int.b - int.a) / (bb - aa);
IF (absError < 100 * NbrRe.Epsilon * absApprox) & (absError > tolerance) THEN
res := RoundoffError
END;
IF NbrRe.Abs( int.b ) < ((1 + 100*NbrRe.Epsilon) * (NbrRe.Abs( int.a ) + 1000/NbrRe.MaxNbr)) THEN
res := RoughIntegrand
END
END Update;
BEGIN
IF f # NIL THEN
res := OKay;
IF a = b THEN error := 0; sign := 1; sumResult := 0
ELSE
maxTol := 0.1; tolerance := NbrRe.Max( 100 * NbrRe.Epsilon, NbrRe.Min( maxTol, error ) );
IF a < b THEN sign := 1; aa := a; bb := b ELSE sign := -1; aa := b; bb := a END;
key := GetKey( aa, bb, bb ); interval := Create( key ); interval.a := aa; interval.b := bb;
Update( interval ); sumResult := interval.soln; sumError := interval.error;
IF NbrCplx.Abs( sumResult ) < 1 THEN error := sumError
ELSE error := NbrCplx.Abs( sumError / sumResult )
END;
IF error > tolerance THEN
NEW( history ); history.Insert( interval, ignor );
LOOP
history.rider.Home; datum := history.rider.Inspect(); interval := datum( CplxInterval );
maxError := interval.error; interval.GetKey( maxKey );
WHILE ~history.rider.eol DO
history.rider.Next; datum := history.rider.Inspect(); interval := datum( CplxInterval );
IF interval.error > maxError THEN
maxError := interval.error; interval.GetKey( maxKey )
END
END;
datum := history.rider.Retrieve( maxKey ); interval := datum( CplxInterval );
history.Delete( maxKey, ignor ); midPoint := (interval.a + interval.b) / 2;
key := GetKey( interval.a, interval.b, midPoint ); intervalL := Create( key );
intervalL.a := interval.a; intervalL.b := midPoint; Update( intervalL );
history.Insert( intervalL, successful );
IF successful & (key # maxKey) THEN
intervalR := Create( maxKey ); intervalR.a := midPoint; intervalR.b := interval.b;
Update( intervalR ); history.Insert( intervalR, ignor )
ELSE
IF successful THEN history.Delete( key, ignor ) END;
history.Insert( interval, ignor ); res := MaxSubDivReached
END;
IF res # OKay THEN EXIT END;
sumResult := sumResult + intervalL.soln + intervalR.soln - interval.soln;
sumError := sumError + intervalL.error + intervalR.error - interval.error;
IF NbrCplx.Abs( sumResult ) < 1 THEN error := sumError
ELSE error := NbrCplx.Abs( sumError / sumResult )
END;
IF error < tolerance THEN EXIT END
END
END
END
ELSE sign := 1; sumResult := 0
END;
RETURN sign * sumResult
END SolveCplx;
PROCEDURE Quadrature;
BEGIN
NEW( node8, 8 ); NEW( wgtGauss4, 4 ); NEW( wgtKronrod8, 8 ); node8[0] := 0.99145537112081263920E0;
node8[1] := 0.949107912; node8[2] := 0.864864423; node8[3] := 0.741531185; node8[4] := 0.586087235;
node8[5] := 0.405845151; node8[6] := 0.207784955; node8[7] := 0.0;
wgtGauss4[0] := 0.129484966; wgtGauss4[1] := 0.279705391; wgtGauss4[2] := 0.381830051;
wgtGauss4[3] := 0.417959183;
wgtKronrod8[0] := 0.022935322;
wgtKronrod8[1] := 0.063092093; wgtKronrod8[2] := 0.104790010; wgtKronrod8[3] := 0.140653260;
wgtKronrod8[4] := 0.169004727; wgtKronrod8[5] := 0.190350578; wgtKronrod8[6] := 0.204432940;
wgtKronrod8[7] := 0.209482141;
NEW( node16, 16 ); NEW( wgtGauss8, 8 ); NEW( wgtKronrod16, 16 );
node16[0] := 0.998002299; node16[1] := 0.987992518; node16[2] := 0.967739076; node16[3] := 0.937273392;
node16[4] := 0.897264532; node16[5] := 0.848206583; node16[6] := 0.790418501; node16[7] := 0.724417731;
node16[8] := 0.650996741; node16[9] := 0.570972173; node16[10] := 0.48508186; node16[11] := 0.39415135;
node16[12] := 0.29918001; node16[13] := 0.20119409; node16[14] := 0.10114207; node16[15] := 0.0;
wgtGauss8[0] := 0.030753242; wgtGauss8[1] := 0.070366047; wgtGauss8[2] := 0.107159220;
wgtGauss8[3] := 0.139570678; wgtGauss8[4] := 0.166269206; wgtGauss8[5] := 0.186161000;
wgtGauss8[6] := 0.198431485; wgtGauss8[7] := 0.202578242;
wgtKronrod16[0] := 0.005377480; wgtKronrod16[1] := 0.015007947; wgtKronrod16[2] := 0.025460847;
wgtKronrod16[3] := 0.035346361; wgtKronrod16[4] := 0.044589751; wgtKronrod16[5] := 0.053481525;
wgtKronrod16[6] := 0.062009568; wgtKronrod16[7] := 0.069854121; wgtKronrod16[8] := 0.076849681;
wgtKronrod16[9] := 0.083080503; wgtKronrod16[10] := 0.088564443; wgtKronrod16[11] := 0.093126598;
wgtKronrod16[12] := 0.096642727; wgtKronrod16[13] := 0.099173599; wgtKronrod16[14] := 0.100769846;
wgtKronrod16[15] := 0.101330007;
NEW( node31, 31 ); NEW( wgtGauss15, 15 ); NEW( wgtKronrod31, 31 );
node31[0] := 0.999484410; node31[1] := 0.996893484; node31[2] := 0.991630997; node31[3] := 0.983668123;
node31[4] := 0.973116323; node31[5] := 0.960021865; node31[6] := 0.944374445; node31[7] := 0.926200047;
node31[8] := 0.905573308; node31[9] := 0.882560536; node31[10] := 0.857205234; node31[11] := 0.829565762;
node31[12] := 0.799727836; node31[13] := 0.767777432; node31[14] := 0.733790062; node31[15] := 0.697850495;
node31[16] := 0.660061064; node31[17] := 0.620526183; node31[18] := 0.579345236; node31[19] := 0.536624148;
node31[20] := 0.492480468; node31[21] := 0.447033770; node31[22] := 0.400401255; node31[23] := 0.352704726;
node31[24] := 0.304073202; node31[25] := 0.254636926; node31[26] := 0.204525117; node31[27] := 0.153869914;
node31[28] := 0.102806938; node31[29] := 0.051471843; node31[30] := 0.0;
wgtGauss15[0] := 0.007968192; wgtGauss15[1] := 0.018466468; wgtGauss15[2] := 0.028784708;
wgtGauss15[3] := 0.038799193; wgtGauss15[4] := 0.048402673; wgtGauss15[5] := 0.057493156;
wgtGauss15[6] := 0.065974230; wgtGauss15[7] := 0.073755975; wgtGauss15[8] := 0.080755895;
wgtGauss15[9] := 0.086899787; wgtGauss15[10] := 0.092122522; wgtGauss15[11] := 0.096368737;
wgtGauss15[12] := 0.099593421; wgtGauss15[13] := 0.101762390; wgtGauss15[14] := 0.102852653;
wgtKronrod31[0] := 0.001389014; wgtKronrod31[1] := 0.003890461; wgtKronrod31[2] := 0.006630704;
wgtKronrod31[3] := 0.009273280; wgtKronrod31[4] := 0.011823015; wgtKronrod31[5] := 0.014369730;
wgtKronrod31[6] := 0.016920889; wgtKronrod31[7] := 0.019414141; wgtKronrod31[8] := 0.021828036;
wgtKronrod31[9] := 0.024191162; wgtKronrod31[10] := 0.026509955; wgtKronrod31[11] := 0.028754049;
wgtKronrod31[12] := 0.030907258; wgtKronrod31[13] := 0.032981447; wgtKronrod31[14] := 0.034979338;
wgtKronrod31[15] := 0.036882365; wgtKronrod31[16] := 0.038678946; wgtKronrod31[17] := 0.040374539;
wgtKronrod31[18] := 0.041969810; wgtKronrod31[19] := 0.043452540; wgtKronrod31[20] := 0.044814800;
wgtKronrod31[21] := 0.046059238; wgtKronrod31[22] := 0.047185547; wgtKronrod31[23] := 0.048185862;
wgtKronrod31[24] := 0.049055435; wgtKronrod31[25] := 0.049795683; wgtKronrod31[26] := 0.050405921;
wgtKronrod31[27] := 0.050881796; wgtKronrod31[28] := 0.051221548; wgtKronrod31[29] := 0.051426129;
wgtKronrod31[30] := 0.051494729
END Quadrature;
BEGIN
Quadrature; MaxIntervals := NbrInt.MaxNbr
END CalcGauss.