175 lines
4.3 KiB
Plaintext
175 lines
4.3 KiB
Plaintext
MODULE numericalIntegrationModula2;
|
|
|
|
(* ISO Modula-2 libraries. *)
|
|
IMPORT LongMath, SLongIO, STextIO;
|
|
|
|
TYPE functionRealToReal = PROCEDURE (LONGREAL) : LONGREAL;
|
|
|
|
PROCEDURE leftRule (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER) : LONGREAL;
|
|
VAR sum : LONGREAL;
|
|
h : LONGREAL;
|
|
i : INTEGER;
|
|
BEGIN
|
|
sum := 0.0;
|
|
h := (b - a) / LFLOAT (n);
|
|
FOR i := 1 TO n DO
|
|
sum := sum + f (a + (h * LFLOAT (i - 1)))
|
|
END;
|
|
RETURN (sum * h)
|
|
END leftRule;
|
|
|
|
PROCEDURE rightRule (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER) : LONGREAL;
|
|
VAR sum : LONGREAL;
|
|
h : LONGREAL;
|
|
i : INTEGER;
|
|
BEGIN
|
|
sum := 0.0;
|
|
h := (b - a) / LFLOAT (n);
|
|
FOR i := 1 TO n DO
|
|
sum := sum + f (a + (h * LFLOAT (i)))
|
|
END;
|
|
RETURN (sum * h)
|
|
END rightRule;
|
|
|
|
PROCEDURE midpointRule (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER) : LONGREAL;
|
|
VAR sum : LONGREAL;
|
|
h : LONGREAL;
|
|
half_h : LONGREAL;
|
|
i : INTEGER;
|
|
BEGIN
|
|
sum := 0.0;
|
|
h := (b - a) / LFLOAT (n);
|
|
half_h := 0.5 * h;
|
|
FOR i := 1 TO n DO
|
|
sum := sum + f (a + (h * LFLOAT (i)) - half_h)
|
|
END;
|
|
RETURN (sum * h)
|
|
END midpointRule;
|
|
|
|
PROCEDURE trapeziumRule (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER) : LONGREAL;
|
|
VAR sum : LONGREAL;
|
|
y0 : LONGREAL;
|
|
y1 : LONGREAL;
|
|
h : LONGREAL;
|
|
i : INTEGER;
|
|
BEGIN
|
|
sum := 0.0;
|
|
h := (b - a) / LFLOAT (n);
|
|
y0 := f (a);
|
|
FOR i := 1 TO n DO
|
|
y1 := f (a + (h * LFLOAT (i)));
|
|
sum := sum + 0.5 * (y0 + y1);
|
|
y0 := y1
|
|
END;
|
|
RETURN (sum * h)
|
|
END trapeziumRule;
|
|
|
|
|
|
PROCEDURE simpsonRule (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER) : LONGREAL;
|
|
VAR sum1 : LONGREAL;
|
|
sum2 : LONGREAL;
|
|
h : LONGREAL;
|
|
half_h : LONGREAL;
|
|
x : LONGREAL;
|
|
i : INTEGER;
|
|
BEGIN
|
|
h := (b - a) / LFLOAT (n);
|
|
half_h := 0.5 * h;
|
|
sum1 := f (a + half_h);
|
|
sum2 := 0.0;
|
|
FOR i := 2 TO n DO
|
|
x := a + (h * LFLOAT (i - 1));
|
|
sum1 := sum1 + f (x + half_h);
|
|
sum2 := sum2 + f (x);
|
|
END;
|
|
RETURN (h / 6.0) * (f (a) + f (b) + (4.0 * sum1) + (2.0 * sum2));
|
|
END simpsonRule;
|
|
|
|
PROCEDURE cube (x : LONGREAL) : LONGREAL;
|
|
BEGIN
|
|
RETURN x * x * x;
|
|
END cube;
|
|
|
|
PROCEDURE reciprocal (x : LONGREAL) : LONGREAL;
|
|
BEGIN
|
|
RETURN 1.0 / x;
|
|
END reciprocal;
|
|
|
|
PROCEDURE identity (x : LONGREAL) : LONGREAL;
|
|
BEGIN
|
|
RETURN x;
|
|
END identity;
|
|
|
|
PROCEDURE printResults (f : functionRealToReal;
|
|
a : LONGREAL;
|
|
b : LONGREAL;
|
|
n : INTEGER;
|
|
nominal : LONGREAL);
|
|
PROCEDURE printOneResult (y : LONGREAL);
|
|
BEGIN
|
|
SLongIO.WriteFloat (y, 16, 20);
|
|
STextIO.WriteString (' (nominal + ');
|
|
SLongIO.WriteFloat (y - nominal, 6, 0);
|
|
STextIO.WriteString (')');
|
|
STextIO.WriteLn;
|
|
END printOneResult;
|
|
BEGIN
|
|
STextIO.WriteString (' left rule ');
|
|
printOneResult (leftRule (f, a, b, n));
|
|
|
|
STextIO.WriteString (' right rule ');
|
|
printOneResult (rightRule (f, a, b, n));
|
|
|
|
STextIO.WriteString (' midpoint rule ');
|
|
printOneResult (midpointRule (f, a, b, n));
|
|
|
|
STextIO.WriteString (' trapezium rule ');
|
|
printOneResult (trapeziumRule (f, a, b, n));
|
|
|
|
STextIO.WriteString (' Simpson rule ');
|
|
printOneResult (simpsonRule (f, a, b, n));
|
|
END printResults;
|
|
|
|
BEGIN
|
|
STextIO.WriteLn;
|
|
|
|
STextIO.WriteString ('x³ in [0,1] with n = 100');
|
|
STextIO.WriteLn;
|
|
printResults (cube, 0.0, 1.0, 100, 0.25);
|
|
|
|
STextIO.WriteLn;
|
|
|
|
STextIO.WriteString ('1/x in [1,100] with n = 1000');
|
|
STextIO.WriteLn;
|
|
printResults (reciprocal, 1.0, 100.0, 1000, LongMath.ln (100.0));
|
|
|
|
STextIO.WriteLn;
|
|
|
|
STextIO.WriteString ('x in [0,5000] with n = 5000000');
|
|
STextIO.WriteLn;
|
|
printResults (identity, 0.0, 5000.0, 5000000, 12500000.0);
|
|
|
|
STextIO.WriteLn;
|
|
|
|
STextIO.WriteString ('x in [0,6000] with n = 6000000');
|
|
STextIO.WriteLn;
|
|
printResults (identity, 0.0, 6000.0, 6000000, 18000000.0);
|
|
|
|
STextIO.WriteLn
|
|
END numericalIntegrationModula2.
|