RosettaCodeData/Task/Numerical-integration/ALGOL-68/numerical-integration.alg

119 lines
3.0 KiB
Plaintext

MODE F = PROC(LONG REAL)LONG REAL;
###############
## left rect ##
###############
PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a;
WHILE x <= b - h DO
sum := sum + (h * f(x));
x +:= h
OD;
sum
END # left rect #;
#################
## right rect ##
#################
PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a + h;
WHILE x <= b DO
sum := sum + (h * f(x));
x +:= h
OD;
sum
END # right rect #;
###############
## mid rect ##
###############
PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a;
WHILE x <= b - h DO
sum := sum + h * f(x + h / 2);
x +:= h
OD;
sum
END # mid rect #;
###############
## trapezium ##
###############
PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= f(a) + f(b);
LONG REAL x:= 1;
WHILE x <= n - 1 DO
sum := sum + 2 * f(a + x * h );
x +:= 1
OD;
(b - a) / (2 * n) * sum
END # trapezium #;
#############
## simpson ##
#############
PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum1:= 0;
LONG REAL sum2:= 0;
INT limit:= n - 1;
FOR i FROM 0 TO limit DO
sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2)
OD;
FOR i FROM 1 TO limit DO
sum2 +:= f(a + h * LONG REAL(i))
OD;
h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END # simpson #;
# test the above procedures #
PROC test integrators = ( STRING legend
, F function
, LONG REAL lower limit
, LONG REAL upper limit
, INT iterations
) VOID:
BEGIN
print( ( legend
, fixed( left rect( function, lower limit, upper limit, iterations ), -20, 6 )
, fixed( right rect( function, lower limit, upper limit, iterations ), -20, 6 )
, fixed( mid rect( function, lower limit, upper limit, iterations ), -20, 6 )
, fixed( trapezium( function, lower limit, upper limit, iterations ), -20, 6 )
, fixed( simpson( function, lower limit, upper limit, iterations ), -20, 6 )
, newline
)
)
END; # test integrators #
print( ( " "
, " left rect"
, " right rect"
, " mid rect"
, " trapezium"
, " simpson"
, newline
)
);
test integrators( "x^3", ( LONG REAL x )LONG REAL: x * x * x, 0, 1, 100 );
test integrators( "1/x", ( LONG REAL x )LONG REAL: 1 / x, 1, 100, 1 000 );
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 5 000, 5 000 000 );
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 6 000, 6 000 000 );
SKIP