Numerical integration: Difference between revisions
→{{header|Nim}}
Line 3,758:
trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</syntaxhighlight>
=={{header|Modula-2}}==
{{works with|GCC|13.1.1}}
For ISO standard Modula-2.
<syntaxhighlight lang="modula2">
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.
</syntaxhighlight>
{{out}}
<pre>$ gm2 -fiso -g -O3 numericalIntegrationModula2.mod && ./a.out
x³ in [0,1] with n = 100
left rule 2.450250000000000E-1 (nominal + -4.97500E-3)
right rule 2.550250000000000E-1 (nominal + 5.02500E-3)
midpoint rule 2.499875000000000E-1 (nominal + -1.25000E-5)
trapezium rule 2.500250000000000E-1 (nominal + 2.50000E-5)
Simpson rule 2.500000000000000E-1 (nominal + -2.71051E-20)
1/x in [1,100] with n = 1000
left rule 4.654991057514676 (nominal + 4.98209E-2)
right rule 4.556981057514676 (nominal + -4.81891E-2)
midpoint rule 4.604762548678375 (nominal + -4.07637E-4)
trapezium rule 4.605986057514676 (nominal + 8.15872E-4)
Simpson rule 4.605170384957142 (nominal + 1.98969E-7)
x in [0,5000] with n = 5000000
left rule 1.249999750000000E+7 (nominal + -2.50000)
right rule 1.250000250000000E+7 (nominal + 2.50000)
midpoint rule 1.250000000000000E+7 (nominal + -1.81899E-12)
trapezium rule 1.250000000000000E+7 (nominal + -1.81899E-12)
Simpson rule 1.250000000000000E+7 (nominal + -9.09495E-13)
x in [0,6000] with n = 6000000
left rule 1.799999700000000E+7 (nominal + -3.00000)
right rule 1.800000300000000E+7 (nominal + 3.00000)
midpoint rule 1.800000000000000E+7 (nominal + 1.81899E-12)
trapezium rule 1.800000000000000E+7 (nominal + 1.81899E-12)
Simpson rule 1.800000000000000E+7 (nominal + 0.00000)
</pre>
=={{header|Nim}}==
| |||