Numerical integration: Difference between revisions
m
→{{header|Wren}}: Minor tidy
(Updated both D entries) |
m (→{{header|Wren}}: Minor tidy) |
||
(112 intermediate revisions by 51 users not shown) | |||
Line 1:
{{task|Arithmetic operations}}
Write functions to calculate the definite integral of a function <big><big> {{math|1=''ƒ(x)''}} </big></big> using ''all'' five of the following methods:
:* [[wp:Rectangle_method|rectangular]]
:** left
:** right
:** midpoint
:* [[wp:Trapezoidal_rule|trapezium]]
:* [[wp:Simpson%27s_rule|Simpson's]]
:** composite
Your functions should take in the upper and lower bounds ({{math|''a''}} and {{math|''b''}}), and the number of approximations to make in that range ({{math|''n''}}).
Assume that your example already has a function that gives values for <big> {{math|1=''ƒ(x)''}} </big>.
Simpson's method is defined by the following pseudo-code:
{| class="mw-collapsible mw-collapsed"
|+ Pseudocode: Simpson's method, composite
|-
|
'''procedure''' quad_simpson_composite(f, a, b, n)
h := (b - a) / n
sum1 := f(a + h/2)
sum2 := 0
loop on i from 1 to (n - 1)
sum1 := sum1 + f(a + h * i + h/2)
sum2 := sum2 + f(a + h * i)
''answer'' := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)
|}
Demonstrate your function by showing the results for:
*
*
*
*
<br/>
;See also:
* [[Active object]] for integrating a function of real time.
* [[Special:PrefixIndex/Numerical integration]] for other integration methods.
<br/>
=={{header|11l}}==
{{trans|Nim}}
<syntaxhighlight lang="11l">F left_rect((Float -> Float) f, Float x, Float h) -> Float
R f(x)
F mid_rect((Float -> Float) f, Float x, Float h) -> Float
R f(x + h / 2)
F right_rect((Float -> Float) f, Float x, Float h) -> Float
R f(x + h)
F trapezium((Float -> Float) f, Float x, Float h) -> Float
R (f(x) + f(x + h)) / 2.0
F simpson((Float -> Float) f, Float x, Float h) -> Float
R (f(x) + 4 * f(x + h / 2) + f(x + h)) / 6.0
F cube(Float x) -> Float
R x * x * x
F reciprocal(Float x) -> Float
R 1 / x
F identity(Float x) -> Float
R x
F integrate(f, a, b, steps, meth)
V h = (b - a) / steps
V ival = h * sum((0 .< steps).map(i -> @meth(@f, @a + i * @h, @h)))
R ival
L(a, b, steps, func, func_name) [(0.0, 1.0, 100, cube, ‘cube’),
(1.0, 100.0, 1000, reciprocal, ‘reciprocal’),
(0.0, 5000.0, 5'000'000, identity, ‘identity’),
(0.0, 6000.0, 6'000'000, identity, ‘identity’)]
L(rule, rule_name) [(left_rect, ‘left_rect’),
(mid_rect, ‘mid_rect’),
(right_rect, ‘right_rect’),
(trapezium, ‘trapezium’),
(simpson, ‘simpson’)]
print("#. integrated using #.\n from #. to #. (#. steps) = #.".format(
func_name, rule_name, a, b, steps, integrate(func, a, b, steps, rule)))</syntaxhighlight>
{{out}}
<pre>
cube integrated using left_rect
from 0 to 1 (100 steps) = 0.245025
cube integrated using mid_rect
from 0 to 1 (100 steps) = 0.2499875
cube integrated using right_rect
from 0 to 1 (100 steps) = 0.255025
cube integrated using trapezium
from 0 to 1 (100 steps) = 0.250025
cube integrated using simpson
from 0 to 1 (100 steps) = 0.25
reciprocal integrated using left_rect
from 1 to 100 (1000 steps) = 4.654991058
reciprocal integrated using mid_rect
from 1 to 100 (1000 steps) = 4.604762549
reciprocal integrated using right_rect
from 1 to 100 (1000 steps) = 4.556981058
reciprocal integrated using trapezium
from 1 to 100 (1000 steps) = 4.605986058
reciprocal integrated using simpson
from 1 to 100 (1000 steps) = 4.605170385
identity integrated using left_rect
from 0 to 5000 (5000000 steps) = 12499997.5
identity integrated using mid_rect
from 0 to 5000 (5000000 steps) = 12500000
identity integrated using right_rect
from 0 to 5000 (5000000 steps) = 12500002.5
identity integrated using trapezium
from 0 to 5000 (5000000 steps) = 12500000
identity integrated using simpson
from 0 to 5000 (5000000 steps) = 12500000
identity integrated using left_rect
from 0 to 6000 (6000000 steps) = 17999997.000000003
identity integrated using mid_rect
from 0 to 6000 (6000000 steps) = 17999999.999999992
identity integrated using right_rect
from 0 to 6000 (6000000 steps) = 18000003.000000003
identity integrated using trapezium
from 0 to 6000 (6000000 steps) = 17999999.999999992
identity integrated using simpson
from 0 to 6000 (6000000 steps) = 17999999.999999992
</pre>
=={{header|ActionScript}}==
Integration functions:
<
{
var sum:Number = 0;
Line 77 ⟶ 185:
}
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
}</
Usage:
<
return (2/(1+ 4*(n*n)));
}
Line 87 ⟶ 195:
trace(trapezium(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
</syntaxhighlight>
=={{header|Ada}}==
Specification of a generic package implementing the five specified kinds of numerical integration:
<
type Scalar is digits <>;
with function F (X : Scalar) return Scalar;
Line 99 ⟶ 208:
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;</
An alternative solution is to pass a function reference to the integration function. This solution is probably slightly faster, and works even with Ada83. One could also make each integration function generic, instead of making the whole package generic.
Body of the package implementing numerical integration:
<
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
H : constant Scalar := (B - A) / Scalar (N);
Line 144 ⟶ 253:
Sum : Scalar := F(A) + F(B);
X : Scalar := 1.0;
begin
while X <= Scalar (N) - 1.0 loop
Sum := Sum + 2.0 * F (A + X * (B - A) / Scalar (N));
Line 154 ⟶ 263:
function Simpsons (A, B : Scalar; N : Positive) return Scalar is
H : constant Scalar := (B - A) / Scalar (N);
begin
for I in
else
Sum_E := Sum_E + F (A + H * Scalar (I));
end if;
end loop;
return (H /
end Simpsons;
end Integrate;</
Test driver:
<
with Integrate;
Line 270 ⟶ 382:
end X;
end Numerical_Integration;
</syntaxhighlight>
=={{header|ALGOL 68}}==
<
###############
Line 357 ⟶ 469:
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</syntaxhighlight>
{{out}}
<pre>
left rect right rect mid rect trapezium simpson
x^3 0.245025 0.255025 0.249988 0.250025 0.250000
1/x 4.654991 4.556981 4.604763 4.605986 4.605170
x 12499997.500000 12500002.500000 12500000.000000 12500000.000000 12500000.000000
x 17999997.000000 18000003.000000 18000000.000000 18000000.000000 18000000.000000</pre>
=={{header|ALGOL W}}==
{{Trans|ALGOL 68}}
<syntaxhighlight lang="algolw">begin % compare some numeric integration methods %
long real procedure leftRect ( long real procedure f
; long real value a, b
; integer value n
) ;
begin
long real h, sum, x;
h := (b - a) / n;
sum := 0;
x := a;
while x <= b - h do begin
sum := sum + (h * f(x));
x := x + h
end;
sum
end leftRect ;
long real procedure rightRect ( long real procedure f
; long real value a, b
; integer value n
) ;
begin
long real h, sum, x;
h := (b - a) / n;
sum := 0;
x := a + h;
while x <= b do begin
sum := sum + (h * f(x));
x := x + h
end;
sum
end rightRect ;
long real procedure midRect ( long real procedure f
; long real value a, b
; integer value n
) ;
begin
long real h, sum, x;
h := (b - a) / n;
sum := 0;
x := a;
while x <= b - h do begin
sum := sum + h * f(x + h / 2);
x := x + h
end;
sum
end midRect ;
long real procedure trapezium ( long real procedure f
; long real value a, b
; integer value n
) ;
begin
long real h, sum, x;
h := (b - a) / n;
sum := f(a) + f(b);
x := 1;
while x <= n - 1 do begin
sum := sum + 2 * f(a + x * h );
x := x + 1
end;
(b - a) / (2 * n) * sum
end trapezium ;
long real procedure simpson ( long real procedure f
; long real value a, b
; integer value n
) ;
begin
long real h, sum1, sum2, x;
integer limit;
h := (b - a) / n;
sum1 := 0;
sum2 := 0;
limit := n - 1;
for i := 0 until limit do sum1 := sum1 + f(a + h * i + h / 2);
for i := 1 until limit do sum2 := sum2 + f(a + h * i);
h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
end simpson ;
% tests the above procedures %
procedure testIntegrators1 ( string(3) value legend
; long real procedure f
; long real value lowerLimit
; long real value upperLimit
; integer value iterations
) ;
write( r_format := "A", r_w := 20, r_d := 6, s_w := 0,
, legend
, leftRect( f, lowerLimit, upperLimit, iterations )
, rightRect( f, lowerLimit, upperLimit, iterations )
, midRect( f, lowerLimit, upperLimit, iterations )
, trapezium( f, lowerLimit, upperLimit, iterations )
, simpson( f, lowerLimit, upperLimit, iterations )
);
procedure testIntegrators2 ( string(3) value legend
; long real procedure f
; long real value lowerLimit
; long real value upperLimit
; integer value iterations
) ;
write( r_format := "A", r_w := 16, r_d := 2, s_w := 0,
, legend
, leftRect( f, lowerLimit, upperLimit, iterations ), " "
, rightRect( f, lowerLimit, upperLimit, iterations ), " "
, midRect( f, lowerLimit, upperLimit, iterations ), " "
, trapezium( f, lowerLimit, upperLimit, iterations ), " "
, simpson( f, lowerLimit, upperLimit, iterations ), " "
);
begin % task test cases %
long real procedure xCubed ( long real value x ) ; x * x * x;
long real procedure oneOverX ( long real value x ) ; 1 / x;
long real procedure xValue ( long real value x ) ; x;
write( " "
, " left rect"
, " right rect"
, " mid rect"
, " trapezium"
, " simpson"
);
testIntegrators1( "x^3", xCubed, 0, 1, 100 );
testIntegrators1( "1/x", oneOverX, 1, 100, 1000 );
testIntegrators2( "x ", xValue, 0, 5000, 5000000 );
testIntegrators2( "x ", xValue, 0, 6000, 6000000 )
end
end.</syntaxhighlight>
=={{header|ATS}}==
<syntaxhighlight lang="ats">
#include "share/atspre_staload.hats"
%{^
#include <math.h>
%}
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
(* This type declarations is for composite quadrature functions for
all the different g0float typekinds. The function must either prove
termination or mask the requirement. (All of ours will prove
termination.) The function to be integrated will not be passed as
an argument, but inlined via the template mechanism. (This design
is more general. It can easily be used to write a quadrature
function that takes the argument, but also can be used for faster
code that requires no function call.) *)
typedef composite_quadrature (tk : tkind) =
(g0float tk, g0float tk, intGte 2) -<> g0float tk
extern fn {tk : tkind}
composite_quadrature$func : g0float tk -<> g0float tk
extern fn {tk : tkind} left_rule : composite_quadrature tk
extern fn {tk : tkind} right_rule : composite_quadrature tk
extern fn {tk : tkind} midpoint_rule : composite_quadrature tk
extern fn {tk : tkind} trapezium_rule : composite_quadrature tk
extern fn {tk : tkind} simpson_rule : composite_quadrature tk
extern fn {tk : tkind}
_one_point_rule$init_x :
g0float tk -<> g0float tk
fn {tk : tkind}
_one_point_rule : composite_quadrature tk =
lam (a, b, n) =>
let
prval [n : int] EQINT () = eqint_make_gint n
macdef f = composite_quadrature$func
val h = (b - a) / g0i2f n
val x0 = _one_point_rule$init_x<tk> h
fun
loop {i : nat | i <= n} .<n - i>.
(i : int i,
sum : g0float tk) :<> g0float tk =
if i = n then
sum
else
loop (succ i, sum + f(x0 + (g0i2f i * h)))
in
loop (0, g0i2f 0) * h
end
(* The left rule, for any floating point type. *)
implement {tk}
left_rule (a, b, n) =
let
implement _one_point_rule$init_x<tk> _ = a
in
_one_point_rule<tk> (a, b, n)
end
(* The right rule, for any floating point type. *)
implement {tk}
right_rule (a, b, n) =
let
implement _one_point_rule$init_x<tk> h = a + h
in
_one_point_rule<tk> (a, b, n)
end
(* The midpoint rule, for any floating point type. *)
implement {tk}
midpoint_rule (a, b, n) =
let
implement _one_point_rule$init_x<tk> h = a + (h / g0i2f 2)
in
_one_point_rule<tk> (a, b, n)
end
implement {tk}
trapezium_rule : composite_quadrature tk =
lam (a, b, n) =>
let
prval [n : int] EQINT () = eqint_make_gint n
macdef f = composite_quadrature$func
val h = (b - a) / g0i2f n
fun
loop {i : pos | i <= n} .<n - i>.
(i : int i,
sum : g0float tk) :<> g0float tk =
if i = n then
sum
else
loop (succ i, sum + f(a + (g0i2f i * h)))
val sum = loop (1, g0i2f 0)
in
((f(a) + sum + sum + f(b)) * h) / g0i2f 2
end
(* Simpson’s 1/3 rule, for any floating point type. *)
implement {tk}
simpson_rule : composite_quadrature tk =
lam (a, b, n) =>
let
(* I have noticed that the Simpson rule is a weighted average of
the trapezium and midpoint rules, which themselves evaluate
the function at different points. Therefore, the following
should be efficient and produce good results. *)
val estimate1 = trapezium_rule<tk> (a, b, n)
val estimate2 = midpoint_rule<tk> (a, b, n)
in
(estimate1 + estimate2 + estimate2) / (g0i2f 3)
end
extern fn {tk : tkind}
fprint_result$rule : composite_quadrature tk
extern fn {tk : tkind}
fprint_result (outf : FILEref,
message : string,
a : g0float tk,
b : g0float tk,
n : intGte 2,
nominal : g0float tk) : void
implement
fprint_result<dblknd> (outf, message, a, b, n, nominal) =
let
val integral = fprint_result$rule<dblknd> (a, b, n)
in
fprint! (outf, " ", message, " ");
ignoret ($extfcall (int, "fprintf", FILEref2star outf,
"%18.15le", integral));
fprint! (outf, " (nominal + ");
ignoret ($extfcall (int, "fprintf", FILEref2star outf,
"% .6le", integral - nominal));
fprint! (outf, ")\n")
end
fn {tk : tkind}
fprint_rule_results (outf : FILEref,
a : g0float tk,
b : g0float tk,
n : intGte 2,
nominal : g0float tk) : void =
let
implement fprint_result$rule<tk> (a, b, n) = left_rule<tk> (a, b, n)
val () = fprint_result (outf, "left rule ", a, b, n, nominal)
implement fprint_result$rule<tk> (a, b, n) = right_rule<tk> (a, b, n)
val () = fprint_result (outf, "right rule ", a, b, n, nominal)
implement fprint_result$rule<tk> (a, b, n) = midpoint_rule<tk> (a, b, n)
val () = fprint_result (outf, "midpoint rule ", a, b, n, nominal)
implement fprint_result$rule<tk> (a, b, n) = trapezium_rule<tk> (a, b, n)
val () = fprint_result (outf, "trapezium rule ", a, b, n, nominal)
implement fprint_result$rule<tk> (a, b, n) = simpson_rule<tk> (a, b, n)
val () = fprint_result (outf, "Simpson rule ", a, b, n, nominal)
in
end
implement
main () =
let
val outf = stdout_ref
val () = fprint! (outf, "\nx³ in [0,1] with n = 100\n")
implement composite_quadrature$func<dblknd> x = x * x * x
val () = fprint_rule_results<dblknd> (outf, 0.0, 1.0, 100, 0.25)
val () = fprint! (outf, "\n1/x in [1,100] with n = 1000\n")
implement composite_quadrature$func<dblknd> x = g0i2f 1 / x
val () = fprint_rule_results<dblknd> (outf, 1.0, 100.0, 1000,
$extfcall (double, "log", 100.0))
val () = fprint! (outf, "\nx in [0,5000] with n = 5000000\n")
implement composite_quadrature$func<dblknd> x = x
val () = fprint_rule_results<dblknd> (outf, 0.0, 5000.0, 5000000,
12500000.0)
val () = fprint! (outf, "\nx in [0,6000] with n = 6000000\n")
implement composite_quadrature$func<dblknd> x = x
val () = fprint_rule_results<dblknd> (outf, 0.0, 6000.0, 6000000,
18000000.0)
val () = fprint! (outf, "\n")
in
0
end
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -Ofast numerical_integration_task.dats -lm && ./a.out
x³ in [0,1] with n = 100
left rule 2.450250000000000e-01 (nominal + -4.975000e-03)
right rule 2.550250000000000e-01 (nominal + 5.025000e-03)
midpoint rule 2.499875000000000e-01 (nominal + -1.250000e-05)
trapezium rule 2.500250000000000e-01 (nominal + 2.500000e-05)
Simpson rule 2.500000000000000e-01 (nominal + 0.000000e+00)
1/x in [1,100] with n = 1000
left rule 4.654991057514675e+00 (nominal + 4.982087e-02)
right rule 4.556981057514675e+00 (nominal + -4.818913e-02)
midpoint rule 4.604762548678376e+00 (nominal + -4.076373e-04)
trapezium rule 4.605986057514674e+00 (nominal + 8.158715e-04)
Simpson rule 4.605170384957142e+00 (nominal + 1.989691e-07)
x in [0,5000] with n = 5000000
left rule 1.249999750000000e+07 (nominal + -2.500000e+00)
right rule 1.250000250000000e+07 (nominal + 2.500000e+00)
midpoint rule 1.250000000000000e+07 (nominal + 0.000000e+00)
trapezium rule 1.250000000000000e+07 (nominal + -1.862645e-09)
Simpson rule 1.250000000000000e+07 (nominal + 0.000000e+00)
x in [0,6000] with n = 6000000
left rule 1.799999700000000e+07 (nominal + -3.000000e+00)
right rule 1.800000300000000e+07 (nominal + 3.000000e+00)
midpoint rule 1.800000000000000e+07 (nominal + 0.000000e+00)
trapezium rule 1.800000000000000e+07 (nominal + 0.000000e+00)
Simpson rule 1.800000000000000e+07 (nominal + 0.000000e+00)
</pre>
=={{header|AutoHotkey}}==
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
<
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
Line 395 ⟶ 905:
fun(x) { ; linear test function
Return x
}</syntaxhighlight>
=={{header|BASIC}}==
{{works with|QuickBasic|4.5}}
{{trans|Java}}
<syntaxhighlight lang="qbasic">FUNCTION leftRect(a, b, n)
h = (b - a) / n
sum = 0
Line 423 ⟶ 931:
h = (b - a) / n
sum = 0
FOR x = a + h / 2 TO b - h / 2 STEP h
sum = sum +
NEXT x
midRect = sum
Line 444 ⟶ 952:
FOR i = 0 TO n-1
sum1 =
NEXT i
Line 452 ⟶ 960:
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION</
=={{header|BBC BASIC}}==
<syntaxhighlight lang="bbcbasic"> *FLOAT64
@% = 12 : REM Column width
PRINT "Function Range L-Rect R-Rect M-Rect Trapeze Simpson"
FOR func% = 1 TO 4
READ x$, l, h, s%
PRINT x$, ; l " - " ; h, FNlrect(x$, l, h, s%) FNrrect(x$, l, h, s%) ;
PRINT FNmrect(x$, l, h, s%) FNtrapeze(x$, l, h, s%) FNsimpson(x$, l, h, s%)
NEXT
END
DATA "x^3", 0, 1, 100
DATA "1/x", 1, 100, 1000
DATA "x", 0, 5000, 5000000
DATA "x", 0, 6000, 6000000
DEF FNlrect(x$, a, b, n%)
LOCAL i%, d, s, x
d = (b - a) / n%
x = a
FOR i% = 1 TO n%
s += d * EVAL(x$)
x += d
NEXT
= s
DEF FNrrect(x$, a, b, n%)
LOCAL i%, d, s, x
d = (b - a) / n%
x = a
FOR i% = 1 TO n%
x += d
s += d * EVAL(x$)
NEXT
= s
DEF FNmrect(x$, a, b, n%)
LOCAL i%, d, s, x
d = (b - a) / n%
x = a
FOR i% = 1 TO n%
x += d/2
s += d * EVAL(x$)
x += d/2
NEXT
= s
DEF FNtrapeze(x$, a, b, n%)
LOCAL i%, d, f, s, x
d = (b - a) / n%
x = b : f = EVAL(x$)
x = a : s = d * (f + EVAL(x$)) / 2
FOR i% = 1 TO n%-1
x += d
s += d * EVAL(x$)
NEXT
= s
DEF FNsimpson(x$, a, b, n%)
LOCAL i%, d, f, s1, s2, x
d = (b - a) / n%
x = b : f = EVAL(x$)
x = a + d/2 : s1 = EVAL(x$)
FOR i% = 1 TO n%-1
x += d/2
s2 += EVAL(x$)
x += d/2
s1 += EVAL(x$)
NEXT
x = a
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</syntaxhighlight>
'''Output:'''
<pre>
Function Range L-Rect R-Rect M-Rect Trapeze Simpson
x^3 0 - 1 0.245025 0.255025 0.2499875 0.250025 0.25
1/x 1 - 100 4.65499106 4.55698106 4.60476255 4.60598606 4.60517038
x 0 - 5000 12499997.5 12500002.5 12500000 12500000 12500000
x 0 - 6000 17999997 18000003 18000000 18000000 18000000
</pre>
=={{header|C}}==
<
#include <stdlib.h>
#include <math.h>
Line 512 ⟶ 1,101:
return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}</
<
double f3(double x)
{
Line 584 ⟶ 1,173:
printf("\n");
}
}</
=={{header|C sharp|C#}}==
<
using System.Collections.Generic;
using System.Linq;
Line 704 ⟶ 1,294:
TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000);
}
}</
Output:
<syntaxhighlight lang="text">0.2499500025
0.24999999875
0.2500500025
Line 725 ⟶ 1,315:
18000003
18000000
18000000</
=={{header|C++}}==
Due to their similarity, it makes sense to make the integration method a policy.
<
template<typename Method, typename F, typename Float>
double integrate(F f, Float a, Float b, int steps, Method m)
Line 784 ⟶ 1,375:
// sample usage
double f(double x) { return x*x;
// inside a function somewhere:
Line 791 ⟶ 1,382:
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double t = integrate(f, 0.0, 1.0, 10, trapezium());
double s = integrate(f, 0.0, 1.0, 10, simpson());</
=={{header|Chapel}}==
<syntaxhighlight lang="chapel">
proc f1(x:real):real {
return x**3;
}
proc f2(x:real):real {
return 1/x;
}
proc f3(x:real):real {
return x;
}
proc leftRectangleIntegration(a: real, b: real, N: int, f): real{
var h: real = (b - a)/N;
var sum: real = 0.0;
var x_n: real;
for n in 0..N-1 {
x_n = a + n * h;
sum = sum + f(x_n);
}
return h * sum;
}
proc rightRectangleIntegration(a: real, b: real, N: int, f): real{
var h: real = (b - a)/N;
var sum: real = 0.0;
var x_n: real;
for n in 0..N-1 {
x_n = a + (n + 1) * h;
sum = sum + f(x_n);
}
return h * sum;
}
proc midpointRectangleIntegration(a: real, b: real, N: int, f): real{
var h: real = (b - a)/N;
var sum: real = 0.0;
var x_n: real;
for n in 0..N-1 {
x_n = a + (n + 0.5) * h;
sum = sum + f(x_n);
}
return h * sum;
}
proc trapezoidIntegration(a: real(64), b: real(64), N: int(64), f): real{
var h: real(64) = (b - a)/N;
var sum: real(64) = f(a) + f(b);
var x_n: real(64);
for n in 1..N-1 {
x_n = a + n * h;
sum = sum + 2.0 * f(x_n);
}
return (h/2.0) * sum;
}
proc simpsonsIntegration(a: real(64), b: real(64), N: int(64), f): real{
var h: real(64) = (b - a)/N;
var sum: real(64) = f(a) + f(b);
var x_n: real(64);
for n in 1..N-1 by 2 {
x_n = a + n * h;
sum = sum + 4.0 * f(x_n);
}
for n in 2..N-2 by 2 {
x_n = a + n * h;
sum = sum + 2.0 * f(x_n);
}
return (h/3.0) * sum;
}
var exact:real;
var calculated:real;
writeln("f(x) = x**3 with 100 steps from 0 to 1");
exact = 0.25;
calculated = leftRectangleIntegration(a = 0.0, b = 1.0, N = 100, f = f1);
writeln("leftRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = rightRectangleIntegration(a = 0.0, b = 1.0, N = 100, f = f1);
writeln("rightRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = midpointRectangleIntegration(a = 0.0, b = 1.0, N = 100, f = f1);
writeln("midpointRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = trapezoidIntegration(a = 0.0, b = 1.0, N = 100, f = f1);
writeln("trapezoidIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = simpsonsIntegration(a = 0.0, b = 1.0, N = 100, f = f1);
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
writeln("f(x) = 1/x with 1000 steps from 1 to 100");
exact = 4.605170;
calculated = leftRectangleIntegration(a = 1.0, b = 100.0, N = 1000, f = f2);
writeln("leftRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = rightRectangleIntegration(a = 1.0, b = 100.0, N = 1000, f = f2);
writeln("rightRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = midpointRectangleIntegration(a = 1.0, b = 100.0, N = 1000, f = f2);
writeln("midpointRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = trapezoidIntegration(a = 1.0, b = 100.0, N = 1000, f = f2);
writeln("trapezoidIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = simpsonsIntegration(a = 1.0, b = 100.0, N = 1000, f = f2);
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
writeln("f(x) = x with 5000000 steps from 0 to 5000");
exact = 12500000;
calculated = leftRectangleIntegration(a = 0.0, b = 5000.0, N = 5000000, f = f3);
writeln("leftRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = rightRectangleIntegration(a = 0.0, b = 5000.0, N = 5000000, f = f3);
writeln("rightRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = midpointRectangleIntegration(a = 0.0, b = 5000.0, N = 5000000, f = f3);
writeln("midpointRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = trapezoidIntegration(a = 0.0, b = 5000.0, N = 5000000, f = f3);
writeln("trapezoidIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = simpsonsIntegration(a = 0.0, b = 5000.0, N = 5000000, f = f3);
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
writeln("f(x) = x with 6000000 steps from 0 to 6000");
exact = 18000000;
calculated = leftRectangleIntegration(a = 0.0, b = 6000.0, N = 6000000, f = f3);
writeln("leftRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = rightRectangleIntegration(a = 0.0, b = 6000.0, N = 6000000, f = f3);
writeln("rightRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = midpointRectangleIntegration(a = 0.0, b = 6000.0, N = 6000000, f = f3);
writeln("midpointRectangleIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = trapezoidIntegration(a = 0.0, b = 6000.0, N = 6000000, f = f3);
writeln("trapezoidIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
calculated = simpsonsIntegration(a = 0.0, b = 6000.0, N = 6000000, f = f3);
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
</syntaxhighlight>
output
<syntaxhighlight lang="text">
f(x) = x**3 with 100 steps from 0 to 1
leftRectangleIntegration: calculated = 0.245025; exact = 0.25; difference = 0.004975
rightRectangleIntegration: calculated = 0.255025; exact = 0.25; difference = 0.005025
midpointRectangleIntegration: calculated = 0.249988; exact = 0.25; difference = 1.25e-05
trapezoidIntegration: calculated = 0.250025; exact = 0.25; difference = 2.5e-05
simpsonsIntegration: calculated = 0.25; exact = 0.25; difference = 5.55112e-17
f(x) = 1/x with 1000 steps from 1 to 100
leftRectangleIntegration: calculated = 4.65499; exact = 4.60517; difference = 0.0498211
rightRectangleIntegration: calculated = 4.55698; exact = 4.60517; difference = 0.0481889
midpointRectangleIntegration: calculated = 4.60476; exact = 4.60517; difference = 0.000407451
trapezoidIntegration: calculated = 4.60599; exact = 4.60517; difference = 0.000816058
simpsonsIntegration: calculated = 4.60517; exact = 4.60517; difference = 3.31627e-06
f(x) = x with 5000000 steps from 0 to 5000
leftRectangleIntegration: calculated = 1.25e+07; exact = 1.25e+07; difference = 2.5
rightRectangleIntegration: calculated = 1.25e+07; exact = 1.25e+07; difference = 2.5
midpointRectangleIntegration: calculated = 1.25e+07; exact = 1.25e+07; difference = 0.0
trapezoidIntegration: calculated = 1.25e+07; exact = 1.25e+07; difference = 1.86265e-09
simpsonsIntegration: calculated = 1.25e+07; exact = 1.25e+07; difference = 3.72529e-09
f(x) = x with 6000000 steps from 0 to 6000
leftRectangleIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 3.0
rightRectangleIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 3.0
midpointRectangleIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 7.45058e-09
trapezoidIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 3.72529e-09
simpsonsIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 0.0
</syntaxhighlight>
=={{header|CoffeeScript}}==
{{trans|python}}
<
rules =
left_rect: (f, x, h) -> f(x)
Line 829 ⟶ 1,583:
result = integrate func, a, b, steps, rule
console.log rule_name, result
</syntaxhighlight>
output
<syntaxhighlight lang="text">
> coffee numerical_integration.coffee
-- tests for cube with 100 steps from 0 to 1
Line 857 ⟶ 1,611:
trapezium 17999999.999999993
simpson 17999999.999999993
</syntaxhighlight>
=={{header|Comal}}==
{{works with|OpenComal on Linux}}
<syntaxhighlight lang="comal">
1000 PRINT "F(X)";" FROM";" TO";" L-Rect";" M-Rect";" R-Rect ";" Trapez";" Simpson"
1010 fromval:=0
1020 toval:=1
1030 PRINT "X^3 ";
1040 PRINT USING "#####": fromval;
1050 PRINT USING "#####": toval;
1060 PRINT USING "###.#########": numint(f1, "L", fromval, toval, 100);
1070 PRINT USING "###.#########": numint(f1, "R", fromval, toval, 100);
1080 PRINT USING "###.#########": numint(f1, "M", fromval, toval, 100);
1090 PRINT USING "###.#########": numint(f1, "T", fromval, toval, 100);
1100 PRINT USING "###.#########": numint(f1, "S", fromval, toval, 100)
1110 //
1120 fromval:=1
1130 toval:=100
1140 PRINT "1/X ";
1150 PRINT USING "#####": fromval;
1160 PRINT USING "#####": toval;
1170 PRINT USING "###.#########": numint(f2, "L", fromval, toval, 1000);
1180 PRINT USING "###.#########": numint(f2, "R", fromval, toval, 1000);
1190 PRINT USING "###.#########": numint(f2, "M", fromval, toval, 1000);
1200 PRINT USING "###.#########": numint(f2, "T", fromval, toval, 1000);
1210 PRINT USING "###.#########": numint(f2, "S", fromval, toval, 1000)
1220 fromval:=0
1230 toval:=5000
1240 PRINT "X ";
1250 PRINT USING "#####": fromval;
1260 PRINT USING "#####": toval;
1270 PRINT USING "#########.###": numint(f3, "L", fromval, toval, 5000000);
1280 PRINT USING "#########.###": numint(f3, "R", fromval, toval, 5000000);
1290 PRINT USING "#########.###": numint(f3, "M", fromval, toval, 5000000);
1300 PRINT USING "#########.###": numint(f3, "T", fromval, toval, 5000000);
1310 PRINT USING "#########.###": numint(f3, "S", fromval, toval, 5000000)
1320 //
1330 fromval:=0
1340 toval:=6000
1350 PRINT "X ";
1360 PRINT USING "#####": fromval;
1370 PRINT USING "#####": toval;
1380 PRINT USING "#########.###": numint(f3, "L", fromval, toval, 6000000);
1390 PRINT USING "#########.###": numint(f3, "R", fromval, toval, 6000000);
1400 PRINT USING "#########.###": numint(f3, "M", fromval, toval, 6000000);
1410 PRINT USING "#########.###": numint(f3, "T", fromval, toval, 6000000);
1420 PRINT USING "#########.###": numint(f3, "S", fromval, toval, 6000000)
1430 END
1440 //
1450 FUNC numint(FUNC f, type$, lbound, rbound, iters) CLOSED
1460 delta:=(rbound-lbound)/iters
1470 integral:=0
1480 CASE type$ OF
1490 WHEN "L", "T", "S"
1500 actval:=lbound
1510 WHEN "M"
1520 actval:=lbound+delta/2
1530 WHEN "R"
1540 actval:=lbound+delta
1550 OTHERWISE
1560 actval:=lbound
1570 ENDCASE
1580 FOR n:=0 TO iters-1 DO
1590 CASE type$ OF
1600 WHEN "L", "M", "R"
1610 integral:+f(actval+n*delta)*delta
1620 WHEN "T"
1630 integral:+delta*(f(actval+n*delta)+f(actval+(n+1)*delta))/2
1640 WHEN "S"
1650 IF n=0 THEN
1660 sum1:=f(lbound+delta/2)
1670 sum2:=0
1680 ELSE
1690 sum1:+f(actval+n*delta+delta/2)
1700 sum2:+f(actval+n*delta)
1710 ENDIF
1720 OTHERWISE
1730 integral:=0
1740 ENDCASE
1750 ENDFOR
1760 IF type$="S" THEN
1770 RETURN (delta/6)*(f(lbound)+f(rbound)+4*sum1+2*sum2)
1780 ELSE
1790 RETURN integral
1800 ENDIF
1810 ENDFUNC
1820 //
1830 FUNC f1(x) CLOSED
1840 RETURN x^3
1850 ENDFUNC
1860 //
1870 FUNC f2(x) CLOSED
1880 RETURN 1/x
1890 ENDFUNC
1900 //
1910 FUNC f3(x) CLOSED
1920 RETURN x
1930 ENDFUNC
</syntaxhighlight>
{{out}}
<pre>
F(X) FROM TO L-Rect M-Rect R-Rect Trapez Simpson
X^3 0 1 0.245025000 0.255025000 0.249987500 0.250025000 0.250000000
1/X 1 100 4.654991058 4.556981058 4.604762549 4.605986058 4.605170385
X 0 5000 12499997.500 12500002.500 12500000.000 12500000.000 12500000.000
X 0 6000 17999997.000 18000003.000 18000000.000 18000000.000 18000000.000
</pre>
=={{header|Common Lisp}}==
<
(* d (loop for x from a below b by d summing (funcall f x))))
Line 887 ⟶ 1,748:
(funcall f b)
(* 4 sum1)
(* 2 sum2))))))</
=={{header|D}}==
<
template integrate(alias method) {
Line 952 ⟶ 1,813:
writeln();
}
}</
Output:
<pre>rectangular left: 0.202500
Line 979 ⟶ 1,840:
===A faster version===
This version avoids function pointers and delegates, same output:
<
template integrate(alias method) {
Line 1,049 ⟶ 1,910:
writeln();
}
}</
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{Trans|Python}}
<syntaxhighlight lang="delphi">program Numerical_integration;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
type
TFx = TFunc<Double, Double>;
TMethod = TFunc<TFx, Double, Double, Double>;
function RectLeft(f: TFx; x, h: Double): Double;
begin
RectLeft := f(x);
end;
function RectMid(f: TFx; x, h: Double): Double;
begin
RectMid := f(x + h / 2);
end;
function RectRight(f: TFx; x, h: Double): Double;
begin
Result := f(x + h);
end;
function Trapezium(f: TFx; x, h: Double): Double;
begin
Result := (f(x) + f(x + h)) / 2.0;
end;
function Simpson(f: TFx; x, h: Double): Double;
begin
Result := (f(x) + 4 * f(x + h / 2) + f(x + h)) / 6.0;
end;
function Integrate(Method: TMethod; f: TFx; a, b: Double; n: Integer): Double;
var
h: Double;
k: integer;
begin
Result := 0;
h := (b - a) / n;
for k := 0 to n - 1 do
Result := Result + Method(f, a + k * h, h);
Result := Result * h;
end;
function f1(x: Double): Double;
begin
Result := x * x * x;
end;
function f2(x: Double): Double;
begin
Result := 1 / x;
end;
function f3(x: Double): Double;
begin
Result := x;
end;
var
fs: array[0..3] of TFx;
mt: array[0..4] of TMethod;
fsNames: array of string = ['x^3', '1/x', 'x', 'x'];
mtNames: array of string = ['RectLeft', 'RectMid', 'RectRight', 'Trapezium', 'Simpson'];
limits: array of array of Double = [[0, 1, 100], [1, 100, 1000], [0, 5000,
5000000], [0, 6000, 6000000]];
i, j, n: integer;
a, b: double;
begin
fs[0] := f1;
fs[1] := f2;
fs[2] := f3;
fs[3] := f3;
mt[0] := RectLeft;
mt[1] := RectMid;
mt[2] := RectRight;
mt[3] := Trapezium;
mt[4] := Simpson;
for i := 0 to High(fs) do
begin
Writeln('Integrate ' + fsNames[i]);
a := limits[i][0];
b := limits[i][1];
n := Trunc(limits[i][2]);
for j := 0 to High(mt) do
Writeln(Format('%.6f', [Integrate(mt[j], fs[i], a, b, n)]));
end;
readln;
end.</syntaxhighlight>
{{out}}
<pre>Integrate x^3
0,245025
0,249988
0,255025
0,250025
0,250000
Integrate 1/x
4,654991
4,604763
4,556981
4,605986
4,605170
Integrate x
12499997,500000
12500000,000000
12500002,500000
12500000,000000
12500000,000000
Integrate x
17999997,000000
18000000,000000
18000003,000000
18000000,000000
18000000,000000</pre>
=={{header|E}}==
{{trans|Python}}
<
def leftRect(f, x, h) {
Line 1,080 ⟶ 2,066:
def h := (b-a) / steps
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
}</
<
# value: 105.33333333333334
? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
# value: 0.10000000002160479</
=={{header|Elixir}}==
<syntaxhighlight lang="elixir">defmodule Numerical do
@funs ~w(leftrect midrect rightrect trapezium simpson)a
def leftrect(f, left,_right), do: f.(left)
def midrect(f, left, right), do: f.((left+right)/2)
def rightrect(f,_left, right), do: f.(right)
def trapezium(f, left, right), do: (f.(left)+f.(right))/2
def simpson(f, left, right), do: (f.(left) + 4*f.((left+right)/2.0) + f.(right)) / 6.0
def integrate(f, a, b, steps) when is_integer(steps) do
delta = (b - a) / steps
Enum.each(@funs, fn fun ->
total = Enum.reduce(0..steps-1, 0, fn i, acc ->
left = a + delta * i
acc + apply(Numerical, fun, [f, left, left+delta])
end)
:io.format "~10s : ~.6f~n", [fun, total * delta]
end)
end
end
f1 = fn x -> x * x * x end
IO.puts "f(x) = x^3, where x is [0,1], with 100 approximations."
Numerical.integrate(f1, 0, 1, 100)
f2 = fn x -> 1 / x end
IO.puts "\nf(x) = 1/x, where x is [1,100], with 1,000 approximations. "
Numerical.integrate(f2, 1, 100, 1000)
f3 = fn x -> x end
IO.puts "\nf(x) = x, where x is [0,5000], with 5,000,000 approximations."
Numerical.integrate(f3, 0, 5000, 5_000_000)
f4 = fn x -> x end
IO.puts "\nf(x) = x, where x is [0,6000], with 6,000,000 approximations."
Numerical.integrate(f4, 0, 6000, 6_000_000)</syntaxhighlight>
{{out}}
<pre>
f(x) = x^3, where x is [0,1], with 100 approximations.
leftrect : 0.245025
midrect : 0.249988
rightrect : 0.255025
trapezium : 0.250025
simpson : 0.250000
f(x) = 1/x, where x is [1,100], with 1,000 approximations.
leftrect : 4.654991
midrect : 4.604763
rightrect : 4.556981
trapezium : 4.605986
simpson : 4.605170
f(x) = x, where x is [0,5000], with 5,000,000 approximations.
leftrect : 12499997.500000
midrect : 12500000.000000
rightrect : 12500002.500000
trapezium : 12500000.000000
simpson : 12500000.000000
f(x) = x, where x is [0,6000], with 6,000,000 approximations.
leftrect : 17999997.000000
midrect : 18000000.000000
rightrect : 18000003.000000
trapezium : 18000000.000000
simpson : 18000000.000000
</pre>
=={{header|Euphoria}}==
<
atom h, sum
h = (bounds[2]-bounds[1])/n
Line 1,167 ⟶ 2,222:
? int_rightrect({0,10},1000,routine_id("x"))
? int_midrect({0,10},1000,routine_id("x"))
? int_simpson({0,10},1000,routine_id("x"))</
Output:
Line 1,185 ⟶ 2,240:
50
</pre>
=={{header|F Sharp}}==
<syntaxhighlight lang="fsharp">
// integration methods
let left f dx x = f x * dx
let right f dx x = f (x + dx) * dx
let mid f dx x = f (x + dx / 2.0) * dx
let trapez f dx x = (f x + f (x + dx)) * dx / 2.0
let simpson f dx x = (f x + 4.0 * f (x + dx / 2.0) + f (x + dx)) * dx / 6.0
// common integration function
let integrate a b f n method =
let dx = (b - a) / float n
[0..n-1] |> Seq.map (fun i -> a + float i * dx) |> Seq.sumBy (method f dx)
// test cases
let methods = [ left; right; mid; trapez; simpson ]
let cases = [
(fun x -> x * x * x), 0.0, 1.0, 100
(fun x -> 1.0 / x), 1.0, 100.0, 1000
(fun x -> x), 0.0, 5000.0, 5000000
(fun x -> x), 0.0, 6000.0, 6000000
]
// execute and output
Seq.allPairs cases methods
|> Seq.map (fun ((f, a, b, n), method) -> integrate a b f n method)
|> Seq.iter (printfn "%f")
</syntaxhighlight>
=={{header|Factor}}==
<syntaxhighlight lang="factor">
USE: math.functions
IN: scratchpad 0 1 [ 3 ^ ] integrate-simpson .
1/4
IN: scratchpad 1000 num-steps set-global
IN: scratchpad 1.0 100 [ -1 ^ ] integrate-simpson .
4.605173316272971
IN: scratchpad 5000000 num-steps set-global
IN: scratchpad 0 5000 [ ] integrate-simpson .
12500000
IN: scratchpad 6000000 num-steps set-global
IN: scratchpad 0 6000 [ ] integrate-simpson .
18000000</syntaxhighlight>
=={{header|Forth}}==
<
defer method ( fn F: x -- fn[x] )
Line 1,226 ⟶ 2,325:
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732</
=={{header|Fortran}}==
In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:
<
real :: elemf, x
elemf = f(x)
end function elemf</
Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function. Methods are collected into a single function in a module.
<
implicit none
Line 1,310 ⟶ 2,409:
end function integrate
end module Integration</
Usage example:
<
use Integration
use FunctionHolder
Line 1,324 ⟶ 2,423:
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')
end program IntegrationTest</
The FunctionHolder module:
<
implicit none
Line 1,340 ⟶ 2,439:
end function afun
end module FunctionHolder</
=={{header|FreeBASIC}}==
Based on the BASIC entry and the BBC BASIC entry
<syntaxhighlight lang="freebasic">' version 17-09-2015
' compile with: fbc -s console
#Define screen_width 1024
#Define screen_height 256
ScreenRes screen_width, screen_height, 8
Width screen_width\8, screen_height\16
Function f1(x As Double) As Double
Return x^3
End Function
Function f2(x As Double) As Double
Return 1/x
End Function
Function f3(x As Double) As Double
Return x
End Function
Function leftrect(a As Double, b As Double, n As Double, _
ByVal f As Function (ByVal As Double) As Double) As Double
Dim As Double sum, x = a, h = (b - a) / n
For i As UInteger = 1 To n
sum = sum + h * f(x)
x = x + h
Next
leftrect = sum
End Function
Function rightrect(a As Double, b As Double, n As Double, _
ByVal f As Function (ByVal As Double) As Double) As Double
Dim As Double sum, x = a, h = (b - a) / n
For i As UInteger = 1 To n
x = x + h
sum = sum + h * f(x)
Next
rightrect = sum
End Function
Function midrect(a As Double, b As Double, n As Double, _
ByVal f As Function (ByVal As Double) As Double) As Double
Dim As Double sum, h = (b - a) / n, x = a + h / 2
For i As UInteger = 1 To n
sum = sum + h * f(x)
x = x + h
Next
midrect = sum
End Function
Function trap(a As Double, b As Double, n As Double, _
ByVal f As Function (ByVal As Double) As Double) As Double
Dim As Double x = a, h = (b - a) / n
Dim As Double sum = h * (f(a) + f(b)) / 2
For i As UInteger = 1 To n -1
x = x + h
sum = sum + h * f(x)
Next
trap = sum
End Function
Function simpson(a As Double, b As Double, n As Double, _
ByVal f As Function (ByVal As Double) As Double) As Double
Dim As UInteger i
Dim As Double sum1, sum2
Dim As Double h = (b - a) / n
For i = 0 To n -1
sum1 = sum1 + f(a + h * i + h / 2)
Next i
For i = 1 To n -1
sum2 = sum2 + f(a + h * i)
Next i
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
End Function
' ------=< main >=------
Dim As Double y
Dim As String frmt = " ##.##########"
Print
Print "function range steps leftrect midrect " + _
"rightrect trap simpson "
Print "f(x) = x^3 0 - 1 100";
Print Using frmt; leftrect(0, 1, 100, @f1); midrect(0, 1, 100, @f1); _
rightrect(0, 1, 100, @f1); trap(0, 1, 100, @f1); simpson(0, 1, 100, @f1)
Print "f(x) = 1/x 1 - 100 1000";
Print Using frmt; leftrect(1, 100, 1000, @f2); midrect(1, 100, 1000, @f2); _
rightrect(1, 100, 1000, @f2); trap(1, 100, 1000, @f2); _
simpson(1, 100, 1000, @f2)
frmt = " #########.###"
Print "f(x) = x 0 - 5000 5000000";
Print Using frmt; leftrect(0, 5000, 5000000, @f3); midrect(0, 5000, 5000000, @f3); _
rightrect(0, 5000, 5000000, @f3); trap(0, 5000, 5000000, @f3); _
simpson(0, 5000, 5000000, @f3)
Print "f(x) = x 0 - 6000 6000000";
Print Using frmt; leftrect(0, 6000, 6000000, @f3); midrect(0, 6000, 6000000, @f3); _
rightrect(0, 6000, 6000000, @f3); trap(0, 6000, 6000000, @f3); _
simpson(0, 6000, 6000000, @f3)
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End</syntaxhighlight>
{{out}}
<pre>function range steps leftrect midrect rightrect trap simpson
f(x) = x^3 0 - 1 100 0.2450250000 0.2499875000 0.2550250000 0.2500250000 0.2500000000
f(x) = 1/x 1 - 100 1000 4.6549910575 4.6047625487 4.5569810575 4.6059860575 4.6051703850
f(x) = x 0 - 5000 5000000 12499997.501 12500000.001 12500002.501 12500000.001 12500000.000
f(x) = x 0 - 6000 6000000 17999997.001 18000000.001 18000003.001 18000000.001 18000000.000</pre>
=={{header|Go}}==
<
import (
Line 1,366 ⟶ 2,600:
spec{0, 6000, 6e6, 18e6, "x", func(x float64) float64 { return x }},
}
// object for associating a printable function name with an integration method
type method struct {
name string
integrate func(spec) float64
}
// integration methods implemented per task description
var methods = []method{
Line 1,381 ⟶ 2,615:
method{"Simpson's ", simpson},
}
func rectLeft(t spec) float64 {
var a adder
r := t.upper - t.lower
nf := float64(t.n)
x0 := t.lower
for i :=
x1 := t.lower + float64(i+1)*r/nf
// x1-x0 better than r/nf.
// (with r/nf, the represenation error accumulates)
x0 = x1
}
return
}
func rectRight(t spec) float64 {
var a adder
r := t.upper - t.lower
nf := float64(t.n)
x0 := t.lower
for i :=
x1 := t.lower + float64(i+1)*r/nf
x0 = x1
}
return
}
func rectMid(t spec) float64 {
var a adder
r := t.upper - t.lower
nf := float64(t.n)
Line 1,421 ⟶ 2,655:
// as well. we just need one extra point, so we use lower-.5.
x0 := t.lower - .5*r/nf
for i :=
x1 := t.lower + (float64(i)+.5)*r/nf
x0 = x1
}
return
}
func trap(t spec) float64 {
var a adder
r := t.upper - t.lower
nf := float64(t.n)
x0 := t.lower
f0 := t.f(x0)
for i :=
x1 := t.lower + float64(i+1)*r/nf
f1 := t.f(x1)
x0, f0 = x1, f1
}
return
}
func simpson(t spec) float64 {
var a adder
r := t.upper - t.lower
nf := float64(t.n)
Line 1,451 ⟶ 2,685:
// we play a little loose with the values used for dx and dx0.
dx0 := r / nf
x0 := t.lower + dx0
for i := 1; i < t.n; i++ {
Line 1,458 ⟶ 2,692:
xmid := (x0 + x1) * .5
dx := x1 - x0
x0 = x1
}
return
}
func sum(v []float64) float64 {
for _, e :=
a.add(e)
}
return a.total()
}
type adder struct {
}
func (a *adder) total() float64 {
return a.sum +
}
func (a *adder) add(x float64) {
sum := a.sum + x
e := sum -
a.e += a.sum - (sum - e) +
}
func main() {
for _, t := range data {
Line 1,513 ⟶ 2,739:
fmt.Println("")
}
}</
{{out}}
<pre>
Integration from 0 to 1 in 100 parts
Exact result 2.5000000e-01 Error
Line 1,555 ⟶ 2,780:
=={{header|Groovy}}==
Solution:
<
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
}
Line 1,598 ⟶ 2,823:
h/3*((fLeft + fRight).sum() + 4*(fMid.sum()))
}
}</
Test:
Each "nRect" (number of rectangles) value given below is the minimum value that meets the tolerance condition for the given circumstances (function-to-integrate, integral-type and integral-bounds).
<
double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d))
Line 1,642 ⟶ 2,867:
assert ((simpsonsIntegral([0d, Math.PI], 1, cubicPoly) - cpIntegralCalc0ToPI)/ cpIntegralCalc0ToPI).abs() < tolerance**2.75 // 1 in 100 billion
double cpIntegralCalcMinusEToPI = (cubicPolyAntiDeriv(Math.PI) - cubicPolyAntiDeriv(-Math.E))
assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</
Requested Demonstrations:
<
println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }])
println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }])
Line 1,675 ⟶ 2,900:
println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }])
println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }])
println ()</
Output:
Line 1,710 ⟶ 2,935:
Different approach from most of the other examples: First, the function ''f'' might be expensive to calculate, and so it should not be evaluated several times. So, ideally, we want to have positions ''x'' and weights ''w'' for each method and then just calculate the approximation of the integral by
<
Second, let's to generalize all integration methods into one scheme. The methods can all be characterized by the coefficients ''vs'' they use in a particular interval. These will be fractions, and for terseness, we extract the denominator as an extra argument ''v''.
Line 1,716 ⟶ 2,941:
Now there are the closed formulas (which include the endpoints) and the open formulas (which exclude them). Let's do the open formulas first, because then the coefficients don't overlap:
<
integrateOpen v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs) * n
Line 1,722 ⟶ 2,947:
ws = concat $ replicate n vs
c = a + h/2
xs = [c + h * fromIntegral i | i <- [0..m-1]]</
Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries:
<
integrateClosed v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs - 1) * n
Line 1,739 ⟶ 2,964:
inter n [] = x : inter (n-1) xs
inter n [y] = (x+y) : inter (n-1) xs
inter n (y:ys) = y : inter n ys</
And now we can just define
<
intRightRect = integrateClosed 1 [0,1]
intMidRect = integrateOpen 1 [1]
intTrapezium = integrateClosed 2 [1,1]
intSimpson = integrateClosed 3 [1,4,1]</
or, as easily, some additional schemes:
<
intOpen1 = integrateOpen 2 [3,3]
intOpen2 = integrateOpen 3 [8,-4,8]</
Some examples:
Line 1,770 ⟶ 2,995:
The whole program:
<syntaxhighlight lang
:: Fractional a
=> (a1 -> a) -> [a1] -> [a] -> a
approx f xs ws =
sum
[ w * f x
| (x, w) <- zip xs ws ]
integrateOpen
:: Fractional a
=> a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateOpen v vs f a b n = approx f xs ws * h / v
where
m = fromIntegral (length vs) * n
ws = concat $ replicate n vs
c = a + h / 2
xs =
[ c + h * fromIntegral i
| i <- [0 .. m - 1] ]
integrateClosed
:: Fractional a
=> a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateClosed v vs f a b n = approx f xs ws * h / v
where
m = fromIntegral (length vs - 1) * n
h = (b - a) / fromIntegral m
[ a + h * fromIntegral i
overlap
:: Num a
=> Int -> [a] -> [a]
overlap n [] = []
overlap n (x:xs) = x : inter n xs
where
inter 1 ys = ys
inter n [] = x : inter (n - 1) xs
inter n [y] = (x + y) : inter (n - 1) xs
inter n (y:ys) = y : inter n ys
uncurry4 :: (t1 -> t2 -> t3 -> t4 -> t) -> (t1, t2, t3, t4) -> t
uncurry4 f ~(a, b, c, d) = f a b c d
-- TEST ----------------------------------------------------------------------
ms
:: Fractional a
=> [(String, (a -> a) -> a -> a -> Int -> a)]
ms =
[ ("rectangular left", integrateClosed 1 [1, 0])
, ("rectangular middle", integrateOpen 1 [1])
, ("rectangular right", integrateClosed 1 [0, 1])
, ("trapezium", integrateClosed 2 [1, 1])
, ("simpson", integrateClosed 3 [1, 4, 1])
]
integrations
:: (Fractional a, Num t, Num t1, Num t2)
=> [(String, (a -> a, t, t1, t2))]
integrations =
[ ("x^3", ((^ 3), 0, 1, 100))
, ("1/x", ((1 /), 1, 100, 1000))
, ("x", (id, 0, 5000, 500000))
, ("x", (id, 0, 6000, 600000))
]
main :: IO ()
main =
mapM_
(\(s, e@(_, a, b, n)) -> do
putStrLn
(concat
[ indent 20 ("f(x) = " ++ s)
, show [a, b]
, " ("
, show n
, " approximations)"
])
mapM_
(\(s, integration) ->
putStrLn (indent 20 (s ++ ":") ++ show (uncurry4 integration e)))
ms
putStrLn [])
integrations
where
indent n = take n . (++ replicate n ' ')</syntaxhighlight>
{{Out}}
<pre>f(x) = x^3 [0.0,1.0] (100 approximations)
rectangular left: 0.24502500000000005
rectangular middle: 0.24998750000000006
rectangular right: 0.25502500000000006
trapezium: 0.25002500000000005
simpson: 0.25000000000000006
f(x) = 1/x [1.0,100.0] (1000 approximations)
rectangular left: 4.65499105751468
rectangular middle: 4.604762548678376
rectangular right: 4.55698105751468
trapezium: 4.605986057514681
simpson: 4.605170384957135
f(x) = x [0.0,5000.0] (500000 approximations)
rectangular left: 1.2499975000000006e7
rectangular middle: 1.2499999999999993e7
rectangular right: 1.2500025000000006e7
trapezium: 1.2500000000000006e7
simpson: 1.2499999999999998e7
f(x) = x [0.0,6000.0] (600000 approximations)
rectangular middle:
rectangular right:
trapezium:
simpson:
Runtime: about 7 seconds.
=={{header|J}}==
===Solution:===
<
'a b steps'=. 3{.y,128
size=. (b - a)%steps
size * +/ u
)
Line 1,879 ⟶ 3,126:
trapezium=: adverb def '-: +/ u y'
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</
===Example usage===
====Required Examples====
<
It=: trapezium integrate
Is=: simpson integrate
Line 1,909 ⟶ 3,156:
1.8e7
] Is 0 6000 6e6
1.8e7</
====Older Examples====
Integrate <code>square</code> (<code>*:</code>) from 0 to π in 10 steps using various methods.
<
10.3095869962
*: trapezium integrate 0 1p1 10
10.3871026879
*: simpson integrate 0 1p1 10
10.3354255601</
Integrate <code>sin</code> from 0 to π in 10 steps using various methods.
<
sin rectangle integrate 0 1p1 10
2.00824840791
Line 1,925 ⟶ 3,172:
1.98352353751
sin simpson integrate 0 1p1 10
2.00000678444</
===Aside===
Note that J has a primitive verb <code>p..</code> for integrating polynomials. For example the integral of <math>x^2</math> (which can be described in terms of its coefficients as <code>0 0 1</code>) is:
<
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3</
That is: <math>0x^0 + 0x^1 + 0x^2 + \tfrac{1}{3}x^3</math><br>
So to integrate <math>x^2</math> from 0 to π :
<
10.3354255601</
That said, J also has <code>d.</code> which can [http://www.jsoftware.com/help/dictionary/dddot.htm integrate] suitable functions.
<
10.3354</
=={{header|Java}}==
<
{
Line 2,066 ⟶ 3,313:
}
}
</syntaxhighlight>
=={{header|jq}}==
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq.'''
The five different integration methods are each presented as independent functions
to facilitate reuse.
<syntaxhighlight lang=jq>
def integrate_left($a; $b; $n; f):
(($b - $a) / $n) as $h
| reduce range(0;$n) as $i (0;
($a + $i * $h) as $x
| . + ($x|f) )
| . * $h;
def integrate_mid($a; $b; $n; f):
(($b - $a) / $n) as $h
| reduce range(0;$n) as $i (0;
($a + $i * $h) as $x
| . + (($x + $h/2) | f) )
| . * $h;
def integrate_right($a; $b; $n; f):
(($b - $a) / $n) as $h
| reduce range(1; $n + 1) as $i (0;
($a + $i * $h) as $x
| . + ($x|f) )
| . * $h;
def integrate_trapezium($a; $b; $n; f):
(($b - $a) / $n) as $h
| reduce range(0;$n) as $i (0;
($a + $i * $h) as $x
| . + ( ($x|f) + (($x + $h)|f)) / 2 )
| . * $h;
def integrate_simpson($a; $b; $n; f):
(($b - $a) / $n) as $h
| reduce range(0;$n) as $i (0;
($a + $i * $h) as $x
| . + ((( ($x|f) + 4 * (($x + ($h/2))|f) + (($x + $h)|f)) / 6)) )
| . * $h;
def demo($a; $b; $n; f):
"Left = \(integrate_left($a;$b;$n;f))",
"Mid = \(integrate_mid ($a;$b;$n;f))",
"Right = \(integrate_right($a;$b;$n;f))",
"Trapezium = \(integrate_trapezium($a;$b;$n;f))",
"Simpson = \(integrate_simpson($a;$b;$n;f))",
"" ;
demo(0; 1; 100; .*.*. ),
demo(1; 100; 1000; 1 / . ),
demo(0; 5000; 5000000; . ),
demo(0; 6000; 6000000; . )
</syntaxhighlight>
{{output}}
<pre>
Left = 0.24502500000000005
Mid = 0.24998750000000006
Right = 0.25502500000000006
Trapezium = 0.250025
Simpson = 0.25
Left = 4.65499105751468
Mid = 4.604762548678376
Right = 4.55698105751468
Trapezium = 4.605986057514676
Simpson = 4.605170384957133
Left = 12499997.5
Mid = 12500000
Right = 12500002.5
Trapezium = 12500000
Simpson = 12500000
Left = 17999997.000000004
Mid = 17999999.999999993
Right = 18000003.000000004
Trapezium = 17999999.999999993
Simpson = 17999999.999999993
</pre>
=={{header|Julia}}==
{{works with|Julia|0.6}}
<syntaxhighlight lang="julia">function simpson(f::Function, a::Number, b::Number, n::Integer)
h = (b - a) / n
s = f(a + h / 2)
for i in 1:(n-1)
s += f(a + h * i + h / 2) + f(a + h * i) / 2
end
return h/6 * (f(a) + f(b) + 4*s)
end
rst =
simpson(x -> x ^ 3, 0, 1, 100),
simpson(x -> 1 / x, 1, 100, 1000),
simpson(x -> x, 0, 5000, 5_000_000),
simpson(x -> x, 0, 6000, 6_000_000)
@show rst</syntaxhighlight>
{{out}}
<pre>rst = (0.25000000000000006, 4.605170384957135, 1.25e7, 1.8e7)</pre>
=={{header|Kotlin}}==
<syntaxhighlight lang="scala">// version 1.1.2
typealias Func = (Double) -> Double
fun integrate(a: Double, b: Double, n: Int, f: Func) {
val h = (b - a) / n
val sum = DoubleArray(5)
for (i in 0 until n) {
val x = a + i * h
sum[0] += f(x)
sum[1] += f(x + h / 2.0)
sum[2] += f(x + h)
sum[3] += (f(x) + f(x + h)) / 2.0
sum[4] += (f(x) + 4.0 * f(x + h / 2.0) + f(x + h)) / 6.0
}
val methods = listOf("LeftRect ", "MidRect ", "RightRect", "Trapezium", "Simpson ")
for (i in 0..4) println("${methods[i]} = ${"%f".format(sum[i] * h)}")
println()
}
fun main(args: Array<String>) {
integrate(0.0, 1.0, 100) { it * it * it }
integrate(1.0, 100.0, 1_000) { 1.0 / it }
integrate(0.0, 5000.0, 5_000_000) { it }
integrate(0.0, 6000.0, 6_000_000) { it }
}</syntaxhighlight>
{{out}}
<pre>
LeftRect = 0.245025
MidRect = 0.249988
RightRect = 0.255025
Trapezium = 0.250025
Simpson = 0.250000
LeftRect = 4.654991
MidRect = 4.604763
RightRect = 4.556981
Trapezium = 4.605986
Simpson = 4.605170
LeftRect = 12499997.500000
MidRect = 12500000.000000
RightRect = 12500002.500000
Trapezium = 12500000.000000
Simpson = 12500000.000000
LeftRect = 17999997.000000
MidRect = 18000000.000000
RightRect = 18000003.000000
Trapezium = 18000000.000000
Simpson = 18000000.000000
</pre>
=={{header|Lambdatalk}}==
Following Python's presentation
<syntaxhighlight lang="scheme">
1) FUNCTIONS
{def left_rect {lambda {:f :x :h} {:f :x}}}
-> left_rect
{def mid_rect {lambda {:f :x :h} {:f {+ :x {/ :h 2}}}}}
-> mid_rect
{def right_rect {lambda {:f :x :h} {:f {+ :x :h}}}}
-> right_rect
{def trapezium {lambda {:f :x :h} {/ {+ {:f :x} {:f {+ :x :h}}} 2}}}
-> trapezium
{def simpson
{lambda {:f :x :h}
{/ {+ {:f :x} {* 4 {:f {+ :x {/ :h 2}}}} {:f {+ :x :h}}} 6}}}
-> simpson
{def cube {lambda {:x} {* :x :x :x}}}
-> cube
{def reciprocal {lambda {:x} {/ 1 :x}}}
-> reciprocal
{def identity {lambda {:x} :x}}
-> identity
{def integrate
{lambda {:f :a :b :steps :meth}
{let { {:f :f} {:a :a} {:steps :steps} {:meth :meth}
{:h {/ {- :b :a} :steps}}
} {* :h {+ {S.map {{lambda {:meth :f :a :h :i}
{:meth :f {+ :a {* :i :h}} :h}
} :meth :f :a :h}
{S.serie 1 :steps}} }}}}}
-> integrate
{def methods left_rect mid_rect right_rect trapezium simpson}
-> methods
2) TESTS
We apply the following template
{b ∫*function* from *a* to *b* steps *steps*}
{table
{tr {td exact value:} {td *value*}} // the awaited value
{S.map {lambda {:m}
{tr {td :m}
{td {integrate *function* *a* *b* *steps* :m}} }}
{methods}} }
to the given *functions* from *a* to *b* with *steps*
and we get:
∫x3 from 0 to 100 steps 100 (computed in 13ms)
exact value: 0.25 // 1/4
left_rect 0.25502500000000006
mid_rect 0.26013825000000007
right_rect 0.26532800000000006
trapezium 0.2601765
simpson 0.260151
∫1/x from 1 to 100 steps 1000 (computed in 94ms)
exact value: 4.605170185988092 // log(100)
left_rect 4.55698105751468
mid_rect 4.511421425235764
right_rect 4.467888185754358
trapezium 4.512434621634517
simpson 4.511759157368674
∫x from 0 to 5000 steps 5000000 (computed in ... 560000m)
exact value: 12500000 // 5000*5000/2
left_rect 12500002.5
mid_rect 12500005
right_rect 12500007.5
trapezium 12500005
simpson 12500005
∫x from 0 to 6000 steps 6000 (computed in 420ms) too impatient for 6000000, sorry
exact value: 18000000 // 6000*6000/2
left_rect 18003000
mid_rect 18006000
right_rect 18009000
trapezium 18006000
simpson 18006000
</syntaxhighlight>
=={{header|Liberty BASIC}}==
Running the big loop value would take a VERY long time & seems unnecessary.<
while 1
read x$
Line 2,187 ⟶ 3,690:
end
</syntaxhighlight>
Numerical integration
Line 2,221 ⟶ 3,724:
Trapezium method 18000000.0000000000 diff 0 0 %
Simpson's Rule 18000000.0000000000 diff 0 0 %
=={{header|Logo}}==
<
output invoke :fn :x
end
Line 2,260 ⟶ 3,761:
print integrate "i.mid "fn2 4 -1 2 ; 2.496091
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</
=={{header|Lua}}==
<
local h = (b - a) / n
local x = a
Line 2,334 ⟶ 3,836:
print( int_methods[i]( function(x) return x end, 0, 5000, 5000000 ) )
print( int_methods[i]( function(x) return x end, 0, 6000, 6000000 ) )
end</
=={{header|
<syntaxhighlight lang="mathematica">leftRect[f_, a_Real, b_Real, N_Integer] :=
Module[{sum = 0, dx = (b - a)/N, x = a, n = N} ,
For[n = N, n > 0, n--, x += dx; sum += f[x];];
Return [ sum*dx ]]
rightRect[f_, a_Real, b_Real, N_Integer] :=
Module[{sum = 0, dx = (b - a)/N, x = a + (b - a)/N, n = N} ,
For[n = N, n > 0, n--, x += dx; sum += f[x];];
Return [ sum*dx ]]
midRect[f_, a_Real, b_Real, N_Integer] :=
Module[{sum = 0, dx = (b - a)/N, x = a + (b - a)/(2 N), n = N} ,
For[n = N, n > 0, n--, x += dx; sum += f[x];];
Return [ sum*dx ]]
trapezium[f_, a_Real, b_Real, N_Integer] :=
Module[{sum = f[a], dx = (b - a)/N, x = a, n = N} ,
For[n = 1, n < N, n++, x += dx; sum += 2 f[x];];
sum += f[b];
Return [ 0.5*sum*dx ]]
simpson[f_, a_Real, b_Real, N_Integer] :=
Module[{sum1 = f[a + (b - a)/(2 N)], sum2 = 0, dx = (b - a)/N, x = a, n = N} ,
For[n = 1, n < N, n++, sum1 += f[a + dx*n + dx/2];
sum2 += f[a + dx*n];];
Return [(dx/6)*(f[a] + f[b] + 4*sum1 + 2*sum2)]]</syntaxhighlight>
<pre>f[x_] := x^3
g[x_] := 1/x
h[x_] := x
Compare[t_] := Apply[ #1, t] & /@ {leftRect, rightRect, midRect, trapezium, simpson}
AccountingForm[
Compare /@ {{f, 0., 1., 100}, {g, 1., 100., 1000},
{h, 0., 5000., 5000000}, {h, 0., 6000., 6000000}}]
->
{{0.255025, 0.265328, 0.260138, 0.250025, 0.25},
{4.55698, 4.46789, 4.51142, 4.60599, 4.60517},
{12500003., 12500008., 12500005., 12500000., 12500000.},
{18000003., 18000009., 18000006., 18000000., 18000000.}}</pre>
=={{header|MATLAB}} / {{header|Octave}}==
For all of the examples given, the function that is passed to the method as parameter f is a function handle.
Function for performing left rectangular integration: leftRectIntegration.m
<
format long;
Line 2,348 ⟶ 3,892:
integral = width * sum( f(x(1:n-1)) );
end</
Function for performing right rectangular integration: rightRectIntegration.m
<
format long;
Line 2,358 ⟶ 3,902:
integral = width * sum( f(x(2:n)) );
end</
Function for performing mid-point rectangular integration: midPointRectIntegration.m
<
format long;
Line 2,368 ⟶ 3,912:
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
end</
Function for performing trapezoidal integration: trapezoidalIntegration.m
<
format long;
Line 2,377 ⟶ 3,921:
integral = trapz( x,f(x) );
end</
Simpson's rule for numerical integration is already included in
<
Using anonymous functions
<
ans =
0.886226925452753</
Using predefined functions
Built-in MATLAB function sin(x):
<
ans =
2.000000000000000</
User defined scripts and functions:
fermiDirac.m
<
k = 8.617343e-5; %Boltazmann's Constant in eV/K
answer = 1./( 1+exp( (x)/(k*2000) ) ); %Fermi-Dirac distribution with mu = 0 and T = 2000K
end</
<
ans =
0.999998006023282</
=={{header|
<syntaxhighlight lang="maxima">right_rect(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n do s: s + subst(x = a + i * h, e),
left_rect(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n do s: s + subst(x = a + (i - 1) * h, e),
mid_rect(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n do s: s + subst(x = a + (i - 1/2) * h, e),
s * h)$
trapezium(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n - 1 do s: s + subst(x = a + i * h, e),
((subst(x = a, e) + subst(x = b, e)) / 2 + s) * h)$
simpson(e, x, a, b, n) := block([h: (b - a) / n, s: 0],
for i from 1 thru n do
s: s + subst(x = a + i * h, e) + 2 * subst(x = a + (i - 1/2) * h, e),
(subst(x = a, e) - subst(x = b, e) + 2 * s) * h / 6)$
/* some tests */
2 * log(2) - 1 - %, bfloat;
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}}==
{{trans|Python}}
<syntaxhighlight lang="nim">type Function = proc(x: float): float
type Rule = proc(f: Function; x, h: float): float
proc leftRect(f: Function; x, h: float): float =
f(x)
proc midRect(f: Function; x, h: float): float =
f(x + h/2.0)
proc rightRect(f: Function; x, h: float): float =
f(x + h)
proc trapezium(f: Function; x, h: float): float =
(f(x) + f(x+h)) / 2.0
proc simpson(f: Function, x, h: float): float =
(f(x) + 4.0*f(x+h/2.0) + f(x+h)) / 6.0
proc cube(x: float): float =
x * x * x
proc reciprocal(x: float): float =
1.0 / x
proc identity(x: float): float =
x
proc integrate(f: Function; a, b: float; steps: int; meth: Rule): float =
let h = (b-a) / float(steps)
for i in 0 ..< steps:
result += meth(f, a+float(i)*h, h)
result = h * result
for fName, a, b, steps, fun in items(
[("cube", 0, 1, 100, cube),
("reciprocal", 1, 100, 1000, reciprocal),
("identity", 0, 5000, 5_000_000, identity),
("identity", 0, 6000, 6_000_000, identity)]):
for rName, rule in items({"leftRect": leftRect, "midRect": midRect,
"rightRect": rightRect, "trapezium": trapezium, "simpson": simpson}):
echo fName, " integrated using ", rName
echo " from ", a, " to ", b, " (", steps, " steps) = ",
integrate(fun, float(a), float(b), steps, rule)</syntaxhighlight>
{{out}}
<pre>cube integrated using leftRect
from 0 to 1 (100 steps) = 0.245025
cube integrated using midRect
from 0 to 1 (100 steps) = 0.2499875000000001
cube integrated using rightRect
from 0 to 1 (100 steps) = 0.2550250000000001
cube integrated using trapezium
from 0 to 1 (100 steps) = 0.250025
cube integrated using simpson
from 0 to 1 (100 steps) = 0.25
reciprocal integrated using leftRect
from 1 to 100 (1000 steps) = 4.65499105751468
reciprocal integrated using midRect
from 1 to 100 (1000 steps) = 4.604762548678376
reciprocal integrated using rightRect
from 1 to 100 (1000 steps) = 4.55698105751468
reciprocal integrated using trapezium
from 1 to 100 (1000 steps) = 4.605986057514676
reciprocal integrated using simpson
from 1 to 100 (1000 steps) = 4.605170384957133
identity integrated using leftRect
from 0 to 5000 (5000000 steps) = 12499997.5
identity integrated using midRect
from 0 to 5000 (5000000 steps) = 12500000.0
identity integrated using rightRect
from 0 to 5000 (5000000 steps) = 12500002.5
identity integrated using trapezium
from 0 to 5000 (5000000 steps) = 12500000.0
identity integrated using simpson
from 0 to 5000 (5000000 steps) = 12500000.0
identity integrated using leftRect
from 0 to 6000 (6000000 steps) = 17999997.0
identity integrated using midRect
from 0 to 6000 (6000000 steps) = 17999999.99999999
identity integrated using rightRect
from 0 to 6000 (6000000 steps) = 18000003.0
identity integrated using trapezium
from 0 to 6000 (6000000 steps) = 17999999.99999999
identity integrated using simpson
from 0 to 6000 (6000000 steps) = 17999999.99999999</pre>
=={{header|OCaml}}==
The problem can be described as integrating using each of a set of methods, over a set of functions, so let us just build the solution in this modular way.
First define the integration function:<syntaxhighlight lang="ocaml">let integrate f a b steps meth =
let h = (b -. a) /. float_of_int steps in
let rec helper i s =
if i >= steps then s
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
in
h *. helper 0 0.</syntaxhighlight>Then list the methods:<syntaxhighlight lang="ocaml">let methods = [
( "rect_l", fun f x _ -> f x);
( "rect_m", fun f x h -> f (x +. h /. 2.) );
( "rect_r", fun f x h -> f (x +. h) );
( "trap", fun f x h -> (f x +. f (x +. h)) /. 2. );
( "simp", fun f x h -> (f x +. 4. *. f (x +. h /. 2.) +. f (x +. h)) /. 6. )
]</syntaxhighlight> and functions (with limits and steps)<syntaxhighlight lang="ocaml">let functions = [
( "cubic", (fun x -> x*.x*.x), 0.0, 1.0, 100);
( "recip", (fun x -> 1.0/.x), 1.0, 100.0, 1000);
( "x to 5e3", (fun x -> x), 0.0, 5000.0, 5_000_000);
( "x to 6e3", (fun x -> x), 0.0, 6000.0, 6_000_000)
]</syntaxhighlight>and finally iterate the integration over both lists:<syntaxhighlight lang="ocaml">let () =
List.iter (fun (s,f,lo,hi,n) ->
Printf.printf "Testing function %s:\n" s;
List.iter (fun (name,meth) ->
Printf.printf " method %s gives %.15g\n" name (integrate f lo hi n meth)
) methods
) functions</syntaxhighlight>Giving the output:
<pre>
Testing function cubic:
method rect_l gives 0.245025
method rect_m gives 0.2499875
method rect_r gives 0.255025
method trap gives 0.250025
method simp gives 0.25
Testing function recip:
method rect_l gives 4.65499105751468
method rect_m gives 4.60476254867838
method rect_r gives 4.55698105751468
method trap gives 4.60598605751468
method simp gives 4.60517038495713
Testing function x to 5e3:
method rect_l gives 12499997.5
method rect_m gives 12500000
method rect_r gives 12500002.5
method trap gives 12500000
method simp gives 12500000
Testing function x to 6e3:
method rect_l gives 17999997
method rect_m gives 18000000
method rect_r gives 18000003
method trap gives 18000000
method simp gives 18000000
</pre>
=={{header|PARI/GP}}==
Note also that double exponential integration is available as <code>intnum(x=a,b,f(x))</code> and Romberg integration is available as <code>intnumromb(x=a,b,f(x))</code>.
<
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
Line 2,477 ⟶ 4,373:
test(x->1/x, 1, 100, 1000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 6000, 6000000)</
Results:
Line 2,510 ⟶ 4,406:
=={{header|Pascal}}==
<
begin
RectLeft := f(xl)
Line 2,550 ⟶ 4,446:
end;
integrate := integral
end;</
=={{header|Perl
{{trans|Raku}}
<syntaxhighlight lang="perl">use feature 'say';
sub leftrect {
my($func, $a, $b, $n) = @_;
my $h = ($b - $a) / $n;
my $sum = 0;
for ($_ = $a; $_ < $b; $_ += $h) { $sum += $func->($_) }
$h * $sum
}
sub rightrect
my($func, $a, $b, $n) = @_;
my $h = ($b - $a) / $n;
my $sum = 0;
for ($_ = $a+$h; $_ < $b+$h; $_ += $h) { $sum += $func->($_) }
$h * $sum
}
sub midrect
my($func, $a, $b, $n) = @_;
my $h = ($b - $a) / $n;
my $sum = 0;
for ($_ = $a + $h/2; $_ < $b; $_ += $h) { $sum += $func->($_) }
$h * $sum
}
sub trapez
my($func, $a, $b, $n) = @_;
my $h = ($b - $a) / $n;
my $
for ($_ = $a+$h; $_ < $b; $_ += $h) { $sum += 2 * $func->($_) }
$h/2 * $sum
}
sub simpsons {
my $h = ($b - $a) / $n;
my $h2 = $h/2;
my $sum1 =
my $sum2 = 0;
for ($_ = $a+$h
$sum1 +=
$sum2 +=
}
}
# round where needed, display in a reasonable format
sub sig {
my $rounded;
if ($value < 10) {
$rounded =~ s/(\.\d*[1-9])0+$/$1/;
$rounded =~ s/\.0+$//;
} else {
$rounded = sprintf "%.1f", $value;
$rounded =~ s/\.0+$//;
}
return $rounded;
}
sub integrate {
my($func, $a, $b, $n, $exact) = @_;
my $f = sub { local $_ = shift; eval $func };
my @res;
push @res, "$func\n in [$a..$b] / $n";
push @res, ' exact result: ' . rnd($exact);
push @res, ' rectangle method left: ' . rnd( leftrect($f, $a, $b, $n));
push @res, ' rectangle method right: ' . rnd(rightrect($f, $a, $b, $n));
push @res, ' rectangle method mid: ' . rnd( midrect($f, $a, $b, $n));
push @res, 'composite trapezoidal rule: ' . rnd( trapez($f, $a, $b, $n));
push @res, ' quadratic simpsons rule: ' . rnd( simpsons($f, $a, $b, $n));
@res;
}
say for integrate('$_ ** 3', 0, 1, 100, 0.25); say '';
say for integrate('1 / $_', 1, 100, 1000, log(100)); say '';
say for integrate('$_', 0, 5_000, 5_000_000, 12_500_000); say '';
say for integrate('$_', 0, 6_000, 6_000_000, 18_000_000);</syntaxhighlight>
{{out}}
<pre>$_ ** 3
in [0..1] / 100
exact result: 0.25
rectangle method left: 0.245025
rectangle method right: 0.255025
rectangle method mid: 0.
composite trapezoidal rule: 0.250025
quadratic simpsons rule: 0.25
1 /
in [1..100] / 1000
exact result: 4.
rectangle method left: 4.
rectangle method right: 4.
rectangle method mid: 4.
composite trapezoidal rule: 4.
quadratic simpsons rule: 4.
in [0..5000] /
exact result: 12500000
rectangle method left:
rectangle method right:
rectangle method mid: 12500000
composite trapezoidal rule: 12500000
quadratic simpsons rule: 12500000
in [0..6000] /
exact result: 18000000
rectangle method left:
rectangle method right:
rectangle method mid: 18000000
composite trapezoidal rule: 18000000
quadratic simpsons rule: 18000000</
=={{header|Phix}}==
<!--<syntaxhighlight lang="phix">(phixonline?)-->
<span style="color: #008080;">function</span> <span style="color: #000000;">rect_left</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000080;font-style:italic;">/*h*/</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rect_mid</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rect_right</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">trapezium</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">))/</span><span style="color: #000000;">2</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">simpson</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">))/</span><span style="color: #000000;">6</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cubed</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">recip</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">x</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">ident</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">x</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">m_id</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">f_id</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">steps</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">accum</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">-</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">steps</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">steps</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">accum</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">m_id</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f_id</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">+</span><span style="color: #000000;">h</span><span style="color: #0000FF;">*</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">*</span><span style="color: #000000;">accum</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">N</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%d"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">N</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%12f"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">N</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1000000</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #008000;">'.'</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">trim_tail</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">trim_tail</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"."</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">name</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">steps</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rid</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">name</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" %-5s %6d - %-5d %10d %12s %12s %12s %12s %12s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">name</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rect_left</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rect_mid</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rect_right</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">trapezium</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">smartp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">integrate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">simpson</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">steps</span><span style="color: #0000FF;">))})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #008000;">"x^3"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">100</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cubed</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #008000;">"1/x"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">100</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">recip</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #008000;">"x"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ident</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #008000;">"x"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ident</span><span style="color: #0000FF;">}}</span>
<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson
x^3 0 - 1 100 0.245025 0.249988 0.255025 0.250025 0.25
1/x 1 - 100 1000 4.654991 4.604763 4.556981 4.605986 4.60517
x 0 - 5000 5000000 12499997.5 12500000 12500002.5 12500000 12500000
x 0 - 6000 6000000 17999997 18000000 18000003 18000000 18000000
</pre>
=={{header|PicoLisp}}==
<
(de leftRect (Fun X)
Line 2,728 ⟶ 4,684:
(*/ H Sum 1.0) ) )
(prinl (round (integrate square 3.0 7.0 30 simpson)))</
Output:
<pre>105.333</pre>
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">integrals: procedure options (main); /* 1 September 2019 */
f: procedure (x, function) returns (float(18));
declare x float(18), function fixed binary;
select (function);
when (1) return (x**3);
when (2) return (1/x);
when (3) return (x);
when (4) return (x);
end;
end f;
declare (a, b) fixed decimal (10);
declare (rect_area, trap_area, Simpson) float(18);
declare (d, dx) float(18);
declare (S1, S2) float(18);
declare N fixed decimal (15), function fixed binary;
declare k fixed decimal (7,2);
put (' Rectangle-left Rectangle-mid Rectangle-right' ||
' Trapezoid Simpson');
do function = 1 to 4;
select(function);
when (1) do; N = 100; a = 0; b = 1; end;
when (2) do; N = 1000; a = 1; b = 100; end;
when (3) do; N = 5000000; a = 0; b = 5000; end;
when (4) do; N = 6000000; a = 0; b = 6000; end;
end;
dx = (b-a)/float(N);
/* Rectangle method, left-side */
rect_area = 0;
do d = 0 to N-1;
rect_area = rect_area + dx*f(a + d*dx, function);
end;
put skip edit (rect_area) (E(25, 15));
/* Rectangle method, mid-point */
rect_area = 0;
do d = 0 to N-1;
rect_area = rect_area + dx*f(a + d*dx + dx/2, function);
end;
put edit (rect_area) (E(25, 15));
/* Rectangle method, right-side */
rect_area = 0;
do d = 1 to N;
rect_area = rect_area + dx*f(a + d*dx, function);
end;
put edit (rect_area) (E(25, 15));
/* Trapezoid method */
trap_area = 0;
do d = 0 to N-1;
trap_area = trap_area + dx*(f(a+d*dx, function) + f(a+(d+1)*dx, function))/2;
end;
put edit (trap_area) (X(1), E(25, 15));
/* Simpson's Rule */
S1 = f(a+dx/2, function);
S2 = 0;
do d = 1 to N-1;
S1 = S1 + f(a+d*dx+dx/2, function);
S2 = S2 + f(a+d*dx, function);
end;
Simpson = dx * (f(a, function) + f(b, function) + 4*S1 + 2*S2) / 6;
put edit (Simpson) (X(1), E(25, 15));
end;
end integrals;
</syntaxhighlight>
<pre>
Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson
2.450250000000000E-0001 2.499875000000000E-0001 2.550250000000000E-0001 2.500250000000000E-0001 2.500000000000000E-0001
4.654991057514676E+0000 4.604762548678375E+0000 4.556981057514676E+0000 4.605986057514676E+0000 4.605170384957142E+0000
1.249999750000000E+0007 1.250000000000000E+0007 1.250000250000000E+0007 1.250000000000000E+0007 1.250000000000000E+0007
1.799999700000000E+0007 1.800000000000000E+0007 1.800000300000000E+0007 1.800000000000000E+0007 1.800000000000000E+0007
</pre>
=={{header|PureBasic}}==
<
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Line 2,835 ⟶ 4,872:
Answer$+"Trapezium="+StrD(Trapezium (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Simpson ="+StrD(Simpson (0,6000,6000000,@Test3()))
MessageRequester("Answer should be 18,000,000",Answer$) </
<pre>Left =0.2353220100
Mid =0.2401367513
Line 2,862 ⟶ 4,899:
=={{header|Python}}==
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output.
<
def left_rect(f,x,h):
Line 2,916 ⟶ 4,953:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</
'''Tests'''
<
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
Line 2,941 ⟶ 4,978:
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</
'''Sample test Output'''
Line 3,026 ⟶ 5,063:
A faster Simpson's rule integrator is
<
h = (b-a)/float(steps)
a1 = a+h/2
s1 = sum( f(a1+i*h) for i in range(0,steps))
s2 = sum( f(a+i*h) for i in range(1,steps))
return (h/6.0)*(f(a)+f(b)+4.0*s1+2.0*s2)</
=={{header|R}}==
The integ function defined below uses arbitrary abscissae and weights passed as argument (resp. u and v). It assumes that f can take a vector argument.
<syntaxhighlight lang="rsplus">integ <- function(f, a, b, n, u, v) {
h <- (b - a) / n
s <- 0
for (i in seq(0, n - 1)) {
s <- s + sum(v * f(a + i * h + u * h))
}
s * h
}
c(rect.left = integ(f, a, b, n, 0, 1),
rect.mid = integ(f, a, b, n, 0.5, 1),
trapezoidal = integ(f, a, b, n, c(0, 1), c(0.5, 0.5)),
simpson = integ(f, a, b, n, c(0, 0.5, 1), c(1, 4, 1) / 6))
}
test(\(x) x^3, 0, 1, 100)
# rect.left rect.right rect.mid trapezoidal simpson
# 0.2450250 0.2550250 0.2499875 0.2500250 0.2500000
test(\(x) 1 / x, 1, 100, 1000)
# rect.left rect.right rect.mid trapezoidal simpson
# 4.654991 4.556981 4.604763 4.605986 4.605170
test(\(x) x, 0, 5000, 5e6)
# rect.left rect.right rect.mid trapezoidal simpson
# 12499998 12500003 12500000 12500000 12500000
test(\(x) x, 0, 6000, 6e6)
# rect.left rect.right rect.mid trapezoidal simpson
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07</syntaxhighlight>
=={{header|Racket}}==
<syntaxhighlight lang="racket">
#lang racket
(define (integrate f a b steps meth)
(define h (/ (- b a) steps))
(* h (for/sum ([i steps])
(meth f (+ a (* h i)) h))))
(define (left-rect f x h) (f x))
(define (mid-rect f x h) (f (+ x (/ h 2))))
(define (right-rect f x h)(f (+ x h)))
(define (trapezium f x h) (/ (+ (f x) (f (+ x h))) 2))
(define (simpson f x h) (/ (+ (f x) (* 4 (f (+ x (/ h 2)))) (f (+ x h))) 6))
(define (test f a b s n)
(displayln n)
(for ([meth (list left-rect mid-rect right-rect trapezium simpson)]
[name '( left-rect mid-rect right-rect trapezium simpson)])
(displayln (~a name ":\t" (integrate f a b s meth))))
(newline))
(test (λ(x) (* x x x)) 0. 1. 100 "CUBED")
(test (λ(x) (/ x)) 1. 100. 1000 "RECIPROCAL")
(test (λ(x) x) 0. 5000. 5000000 "IDENTITY")
(test (λ(x) x) 0. 6000. 6000000 "IDENTITY")
</syntaxhighlight>
Output:
<syntaxhighlight lang="racket">
CUBED
left-rect: 0.24502500000000005
mid-rect: 0.24998750000000006
right-rect: 0.25502500000000006
trapezium: 0.250025
simpson: 0.25
RECIPROCAL
left-rect: 4.65499105751468
mid-rect: 4.604762548678376
right-rect: 4.55698105751468
trapezium: 4.605986057514676
simpson: 4.605170384957133
IDENTITY
left-rect: 12499997.5
mid-rect: 12500000.0
right-rect: 12500002.5
trapezium: 12500000.0
simpson: 12500000.0
IDENTITY
left-rect: 17999997.000000004
mid-rect: 17999999.999999993
right-rect: 18000003.000000004
trapezium: 17999999.999999993
simpson: 17999999.999999993
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
The addition of <tt>'''Promise'''</tt>/<tt>'''await'''</tt> allows for concurrent computation, and brings a significant speed-up in running time. Which is not to say that it makes this code fast, but it does make it less slow.
Note that these integrations are done with rationals rather than floats, so should be fairly precise (though of course with so few iterations they are not terribly accurate (except when they are)). Some of the sums do overflow into <tt>Num</tt> (floating point)--currently Rakudo allows 64-bit denominators--but at least all of the interval arithmetic is exact.
{{works with|Rakudo|2018.09}}
<syntaxhighlight lang="raku" line>use MONKEY-SEE-NO-EVAL;
sub leftrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $end = $b-$h;
my $sum = 0;
loop (my $i = $a; $i <= $end; $i += $h) { $sum += f($i) }
$h * $sum;
}
sub rightrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $sum = 0;
loop (my $i = $a+$h; $i <= $b; $i += $h) { $sum += f($i) }
}
sub midrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $sum = 0;
my ($start, $end) = $a+$h/2, $b-$h/2;
loop (my $i = $start; $i <= $end; $i += $h) { $sum += f($i) }
$h * $sum;
}
sub trapez(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $partial-sum = 0;
my ($start, $end) = $a+$h, $b-$h;
loop (my $i = $start; $i <= $end; $i += $h) { $partial-sum += f($i) * 2 }
$h / 2 * ( f($a) + f($b) + $partial-sum );
}
sub simpsons(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $h2 = $h/2;
my ($start, $end) = $a+$h, $b-$h;
my $sum1 = f($a + $h2);
my $sum2 = 0;
loop (my $i = $start; $i <= $end; $i += $h) {
$sum1 += f($i + $h2);
$sum2 += f($i);
}
($h / 6) * (f($a) + f($b) + 4*$sum1 + 2*$sum2);
}
sub integrate($f, $a, $b, $n, $exact) {
my $e = 0.000001;
my $r0 = "$f\n in [$a..$b] / $n\n"
~ ' exact result: '~ $exact.round($e);
my ($r1,$r2,$r3,$r4,$r5);
my &f;
EVAL "&f = $f";
my $p1 = Promise.start( { $r1 = ' rectangle method left: '~ leftrect(&f, $a, $b, $n).round($e) } );
my $p2 = Promise.start( { $r2 = ' rectangle method right: '~ rightrect(&f, $a, $b, $n).round($e) } );
my $p3 = Promise.start( { $r3 = ' rectangle method mid: '~ midrect(&f, $a, $b, $n).round($e) } );
my $p4 = Promise.start( { $r4 = 'composite trapezoidal rule: '~ trapez(&f, $a, $b, $n).round($e) } );
my $p5 = Promise.start( { $r5 = ' quadratic simpsons rule: '~ simpsons(&f, $a, $b, $n).round($e) } );
await $p1, $p2, $p3, $p4, $p5;
$r0, $r1, $r2, $r3, $r4, $r5;
}
.say for integrate '{ $_ ** 3 }', 0, 1, 100, 0.25; say '';
.say for integrate '1 / *', 1, 100, 1000, log(100); say '';
.say for integrate '*.self', 0, 5_000, 5_000_000, 12_500_000; say '';
.say for integrate '*.self', 0, 6_000, 6_000_000, 18_000_000;</syntaxhighlight>
{{out}}
<pre>{ $_ ** 3 }
in [0..1] / 100
exact result: 0.25
rectangle method left: 0.245025
rectangle method right: 0.255025
rectangle method mid: 0.249988
composite trapezoidal rule: 0.250025
quadratic simpsons rule: 0.25
1 / *
in [1..100] / 1000
exact result: 4.60517
rectangle method left: 4.654991
rectangle method right: 4.556981
rectangle method mid: 4.604763
composite trapezoidal rule: 4.605986
quadratic simpsons rule: 4.60517
*.self
in [0..5000] / 5000000
exact result: 12500000
rectangle method left: 12499997.5
rectangle method right: 12500002.5
rectangle method mid: 12500000
composite trapezoidal rule: 12500000
quadratic simpsons rule: 12500000
*.self
in [0..6000] / 6000000
exact result: 18000000
rectangle method left: 17999997
rectangle method right: 18000003
rectangle method mid: 18000000
composite trapezoidal rule: 18000000
quadratic simpsons rule: 18000000</pre>
=={{header|REXX}}==
Note: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''.
<syntaxhighlight lang="rexx">/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
numeric digits 20 /*use twenty decimal digits precision. */
do test=1 for 4; say /*perform the 4 different test suites. */
if test==1 then do; L= 0; H= 1; i= 100; end
if test==2 then do; L= 1; H= 100; i= 1000; end
if test==3 then do; L= 0; H= 5000; i= 5000000; end
if test==4 then do; L= 0; H= 6000; i= 6000000; end
say center('test' test, 79, "═") /*display a header for the test suite. */
say ' left rectangular('L", "H', 'i") ──► " left_rect(L, H, i)
say ' midpoint rectangular('L", "H', 'i") ──► " midpoint_rect(L, H, i)
say ' right rectangular('L", "H', 'i") ──► " right_rect(L, H, i)
say ' Simpson('L", "H', 'i") ──► " Simpson(L, H, i)
say ' trapezium('L", "H', 'i") ──► " trapezium(L, H, i)
end /*test*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
f: parse arg y; if test>2 then return y /*choose the "as─is" function. */
if test==1 then return y**3 /* " " cube function. */
return 1/y /* " " reciprocal " */
/*──────────────────────────────────────────────────────────────────────────────────────*/
left_rect: procedure expose test; parse arg a,b,#; $= 0; h= (b-a)/#
do x=a by h for #; $= $ + f(x)
end /*x*/
return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
midpoint_rect: procedure expose test; parse arg a,b,#; $= 0; h= (b-a)/#
do x=a+h/2 by h for #; $= $ + f(x)
end /*x*/
return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
right_rect: procedure expose test; parse arg a,b,#; $= 0; h= (b-a)/#
do x=a+h by h for #; $= $ + f(x)
end /*x*/
return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
Simpson: procedure expose test; parse arg a,b,#; h= (b-a)/#
hh= h/2; $= f(a + hh)
@= 0; do x=1 for #-1; hx=h*x + a; @= @ + f(hx)
$= $ + f(hx + hh)
end /*x*/
return h * (f(a) + f(b) + 4*$ + 2*@) / 6
/*──────────────────────────────────────────────────────────────────────────────────────*/
trapezium: procedure expose test; parse arg a,b,#; $= 0; h= (b-a)/#
do x=a by h for #; $= $ + (f(x) + f(x+h))
end /*x*/
return $*h/2</syntaxhighlight>
{{out|output|text= when using the default inputs:}}
<pre>
════════════════════════════════════test 1═════════════════════════════════════
left rectangular(0, 1, 100) ──► 0.245025
midpoint rectangular(0, 1, 100) ──► 0.2499875
right rectangular(0, 1, 100) ──► 0.255025
Simpson(0, 1, 100) ──► 0.25
trapezium(0, 1, 100) ──► 0.250025
════════════════════════════════════test 2═════════════════════════════════════
left rectangular(1, 100, 1000) ──► 4.6549910575146761473
midpoint rectangular(1, 100, 1000) ──► 4.604762548678375185
right rectangular(1, 100, 1000) ──► 4.5569810575146761472
Simpson(1, 100, 1000) ──► 4.6051703849571421725
trapezium(1, 100, 1000) ──► 4.605986057514676146
════════════════════════════════════test 3═════════════════════════════════════
left rectangular(0, 5000, 5000000) ──► 12499997.5
midpoint rectangular(0, 5000, 5000000) ──► 12500000
right rectangular(0, 5000, 5000000) ──► 12500002.5
Simpson(0, 5000, 5000000) ──► 12500000
trapezium(0, 5000, 5000000) ──► 12500000
════════════════════════════════════test 4═════════════════════════════════════
left rectangular(0, 6000, 6000000) ──► 17999997
midpoint rectangular(0, 6000, 6000000) ──► 18000000
right rectangular(0, 6000, 6000000) ──► 18000003
Simpson(0, 6000, 6000000) ──► 18000000
trapezium(0, 6000, 6000000) ──► 18000000
</pre>
=={{header|Ring}}==
<syntaxhighlight lang="ring">
# Project : Numerical integration
decimals(8)
data = [["pow(x,3)",0,1,100], ["1/x",1, 100,1000], ["x",0,5000,5000000], ["x",0,6000,6000000]]
see "Function Range L-Rect R-Rect M-Rect Trapeze Simpson" + nl
for p = 1 to 4
d1 = data[p][1]
d2 = data[p][2]
d3 = data[p][3]
d4 = data[p][4]
see "" + d1 + " " + d2 + " - " + d3 + " " + lrect(d1, d2, d3, d4) + " " + rrect(d1, d2, d3, d4)
see " " + mrect(d1, d2, d3, d4) + " " + trapeze(d1, d2, d3, d4) + " " + simpson(d1, d2, d3, d4) + nl
next
func lrect(x2, a, b, n)
s = 0
d = (b - a) / n
x = a
for i = 1 to n
eval("result = " + x2)
s = s + d * result
x = x + d
next
return s
func rrect(x2, a, b, n)
s = 0
d = (b - a) / n
x = a
for i = 1 to n
x = x + d
eval("result = " + x2)
s = s + d *result
next
return s
func mrect(x2, a, b, n)
s = 0
d = (b - a) / n
x = a
for i = 1 to n
x = x + d/2
eval("result = " + x2)
s = s + d * result
x = x +d/2
next
return s
func trapeze(x2, a, b, n)
s = 0
d = (b - a) / n
x = b
eval("result = " + x2)
f = result
x = a
eval("result = " + x2)
s = d * (f + result) / 2
for i = 1 to n-1
x = x + d
eval("result = " + x2)
s = s + d * result
next
return s
func simpson(x2, a, b, n)
s1 = 0
s = 0
d = (b - a) / n
x = b
eval("result = " + x2)
f = result
x = a + d/2
eval("result = " + x2)
s1 = result
for i = 1 to n-1
x = x + d/2
eval("result = " + x2)
s = s + result
x = x + d/2
eval("result = " + x2)
s1 = s1 + result
next
x = a
eval("result = " + x2)
return (d / 6) * (f + result + 4 * s1 + 2 * s)
</syntaxhighlight>
Output:
<pre>
Function Range L-Rect R-Rect M-Rect Trapeze Simpson
pow(x,3) 0 - 1 0.245025 0.255025 0.2499875 0.250025 0.25
1/x 1 - 100 4.65499106 4.55698106 4.60476255 4.60598606 4.60517038
x 0 - 5000 12499997.5 12500002.5 12500000 12500000 12500000
x 0 - 6000 17999997 18000003 18000000 18000000 18000000
</pre>
=={{header|Ruby}}==
{{trans|Tcl}}
<
f.call(left)
end
Line 3,255 ⟶ 5,512:
printf " %-10s %s\t(%.1f%%)\n", method, int, diff
end
end</
outputs
<pre>integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
Line 3,269 ⟶ 5,526:
trapezium 1.98352353750945 (-0.8%)
simpson 2.0000067844418 (0.0%)</pre>
=={{header|Rust}}==
This is a partial solution and only implements trapezium integration.
<syntaxhighlight lang="rust">fn integral<F>(f: F, range: std::ops::Range<f64>, n_steps: u32) -> f64
where F: Fn(f64) -> f64
{
let step_size = (range.end - range.start)/n_steps as f64;
let mut integral = (f(range.start) + f(range.end))/2.;
let mut pos = range.start + step_size;
while pos < range.end {
integral += f(pos);
pos += step_size;
}
integral * step_size
}
fn main() {
println!("{}", integral(|x| x.powi(3), 0.0..1.0, 100));
println!("{}", integral(|x| 1.0/x, 1.0..100.0, 1000));
println!("{}", integral(|x| x, 0.0..5000.0, 5_000_000));
println!("{}", integral(|x| x, 0.0..6000.0, 6_000_000));
}</syntaxhighlight>
{{out}}
<pre>0.2500250000000004
4.605986057514688
12500000.000728702
18000000.001390498</pre>
=={{header|Scala}}==
<
def leftRect(f:Double=>Double, a:Double, b:Double)=f(a)
def midRect(f:Double=>Double, a:Double, b:Double)=f((a+b)/2)
Line 3,305 ⟶ 5,591:
print(fn3, 0, 6000, 6000000)
}
}</
Output:
<pre>rectangular left : 0,245025
Line 3,333 ⟶ 5,619:
=={{header|Scheme}}==
<
(define h (/ (- b a) steps))
(* h
Line 3,353 ⟶ 5,639:
(define rr (integrate square 0 1 10 right-rect))
(define t (integrate square 0 1 10 trapezium))
(define s (integrate square 0 1 10 simpson))</
=={{header|SequenceL}}==
<syntaxhighlight lang="sequencel">import <Utilities/Conversion.sl>;
import <Utilities/Sequence.sl>;
integrateLeft(f, a, b, n) :=
let
h := (b - a) / n;
vals[x] := f(x) foreach x within (0 ... (n-1)) * h + a;
in
h * sum(vals);
integrateRight(f, a, b, n) :=
let
h := (b - a) / n;
vals[x] := f(x+h) foreach x within (0 ... (n-1)) * h + a;
in
h * sum(vals);
integrateMidpoint(f, a, b, n) :=
let
h := (b - a) / n;
vals[x] := f(x+h/2.0) foreach x within (0 ... (n-1)) * h + a;
in
h * sum(vals);
integrateTrapezium(f, a, b, n) :=
let
h := (b - a) / n;
vals[i] := 2.0 * f(a + i * h) foreach i within 1 ... n-1;
in
h * (sum(vals) + f(a) + f(b)) / 2.0;
integrateSimpsons(f, a, b, n) :=
let
h := (b - a) / n;
vals1[i] := f(a + h * i + h / 2.0) foreach i within 0 ... n-1;
vals2[i] := f(a + h * i) foreach i within 1 ... n-1;
in
h / 6.0 * (f(a) + f(b) + 4.0 * sum(vals1) + 2.0 * sum(vals2));
xCubed(x) := x^3;
xInverse(x) := 1/x;
identity(x) := x;
tests[method] :=
[method(xCubed, 0.0, 1.0, 100),
method(xInverse, 1.0, 100.0, 1000),
method(identity, 0.0, 5000.0, 5000000),
method(identity, 0.0, 6000.0, 6000000)]
foreach method within [integrateLeft, integrateRight, integrateMidpoint, integrateTrapezium, integrateSimpsons];
//String manipulation for ouput display.
main :=
let
heading := [["Func", "Range\t", "L-Rect\t", "R-Rect\t", "M-Rect\t", "Trapezium", "Simpson"]];
ranges := [["0 - 1\t", "1 - 100\t", "0 - 5000", "0 - 6000"]];
funcs := [["x^3", "1/x", "x", "x"]];
in
delimit(delimit(heading ++ transpose(funcs ++ ranges ++ trimEndZeroes(floatToString(tests, 8))), '\t'), '\n');
trimEndZeroes(x(1)) := x when size(x) = 0 else x when x[size(x)] /= '0' else trimEndZeroes(x[1...size(x)-1]);</syntaxhighlight>
{{out}}
<pre style="height: 25ex; overflow: scroll">
"Func Range L-Rect R-Rect M-Rect Trapezium Simpson
x^3 0 - 1 0.245025 0.255025 0.2499875 0.250025 0.25
1/x 1 - 100 4.65499106 4.55698106 4.60476255 4.60598606 4.60517038
x 0 - 5000 12499997.5 12500002.5 12500000. 12500000. 12500000.
x 0 - 6000 17999997. 18000003. 18000000. 18000000. 18000000."
</pre>
=={{header|Sidef}}==
{{trans|Raku}}
<syntaxhighlight lang="ruby">func sum(f, start, from, to) {
var s = 0;
RangeNum(start, to, from-start).each { |i|
s += f(i);
}
return s
}
func leftrect(f, a, b, n) {
var h = ((b - a) / n);
h * sum(f, a, a+h, b-h);
}
func rightrect(f, a, b, n) {
var h = ((b - a) / n);
h * sum(f, a+h, a + 2*h, b);
}
func midrect(f, a, b, n) {
var h = ((b - a) / n);
h * sum(f, a + h/2, a + h + h/2, b - h/2)
}
func trapez(f, a, b, n) {
var h = ((b - a) / n);
h/2 * (f(a) + f(b) + sum({ f(_)*2 }, a+h, a + 2*h, b-h));
}
func simpsons(f, a, b, n) {
var h = ((b - a) / n);
var h2 = h/2;
var sum1 = f(a + h2);
var sum2 = 0;
sum({|i| sum1 += f(i + h2); sum2 += f(i); 0 }, a+h, a+h+h, b-h);
h/6 * (f(a) + f(b) + 4*sum1 + 2*sum2);
}
func tryem(label, f, a, b, n, exact) {
say "\n#{label}\n in [#{a}..#{b}] / #{n}";
say(' exact result: ', exact);
say(' rectangle method left: ', leftrect(f, a, b, n));
say(' rectangle method right: ', rightrect(f, a, b, n));
say(' rectangle method mid: ', midrect(f, a, b, n));
say('composite trapezoidal rule: ', trapez(f, a, b, n));
say(' quadratic simpsons rule: ', simpsons(f, a, b, n));
}
tryem('x^3', { _ ** 3 }, 0, 1, 100, 0.25);
tryem('1/x', { 1 / _ }, 1, 100, 1000, log(100));
tryem('x', { _ }, 0, 5_000, 5_000_000, 12_500_000);
tryem('x', { _ }, 0, 6_000, 6_000_000, 18_000_000);</syntaxhighlight>
=={{header|Standard ML}}==
<
val h = (b - a) / real steps
fun helper (i, s) =
Line 3,378 ⟶ 5,792:
val rr = integrate (square, 0.0, 1.0, 10, right_rect)
val t = integrate (square, 0.0, 1.0, 10, trapezium )
val s = integrate (square, 0.0, 1.0, 10, simpson )</
=={{header|Stata}}==
<syntaxhighlight lang="text">mata
function integrate(f,a,b,n,u,v) {
s = 0
h = (b-a)/n
m = length(u)
for (i=0; i<n; i++) {
x = a+i*h
for (j=1; j<=m; j++) s = s+v[j]*(*f)(x+h*u[j])
}
return(s*h)
}
function log_(x) {
return(log(x))
}
function id(x) {
return(x)
}
function cube(x) {
return(x*x*x)
}
function inv(x) {
return(1/x)
}
function test(f,a,b,n) {
return(integrate(f,a,b,n,(0,1),(1,0)),
integrate(f,a,b,n,(0,1),(0,1)),
integrate(f,a,b,n,(0.5),(1)),
integrate(f,a,b,n,(0,1),(0.5,0.5)),
integrate(f,a,b,n,(0,1/2,1),(1/6,4/6,1/6)))
}
test(&cube(),0,1,100)
test(&inv(),1,100,1000)
test(&id(),0,5000,5000000)
test(&id(),0,6000,6000000)
end</syntaxhighlight>
'''Output'''
<pre> 1 2 3 4 5
+--------------------------------------------------------+
1 | .245025 .255025 .2499875 .250025 .25 |
+--------------------------------------------------------+
1 2 3 4 5
+-----------------------------------------------------------------------+
1 | 4.654991058 4.556981058 4.604762549 4.605986058 4.605170385 |
+-----------------------------------------------------------------------+
1 2 3 4 5
+------------------------------------------------------------------+
1 | 12499997.5 12500002.5 12500000 12500000 12500000 |
+------------------------------------------------------------------+
1 2 3 4 5
+--------------------------------------------------------+
1 | 17999997 18000003 18000000 18000000 18000000 |
+--------------------------------------------------------+</pre>
=={{header|Swift}}==
<syntaxhighlight lang="swift">public enum IntegrationType : CaseIterable {
case rectangularLeft
case rectangularRight
case rectangularMidpoint
case trapezium
case simpson
}
public func integrate(
from: Double,
to: Double,
n: Int,
using: IntegrationType = .simpson,
f: (Double) -> Double
) -> Double {
let integrationFunc: (Double, Double, Int, (Double) -> Double) -> Double
switch using {
case .rectangularLeft:
integrationFunc = integrateRectL
case .rectangularRight:
integrationFunc = integrateRectR
case .rectangularMidpoint:
integrationFunc = integrateRectMid
case .trapezium:
integrationFunc = integrateTrapezium
case .simpson:
integrationFunc = integrateSimpson
}
return integrationFunc(from, to, n, f)
}
private func integrateRectL(from: Double, to: Double, n: Int, f: (Double) -> Double) -> Double {
let h = (to - from) / Double(n)
var x = from
var sum = 0.0
while x <= to - h {
sum += f(x)
x += h
}
return h * sum
}
private func integrateRectR(from: Double, to: Double, n: Int, f: (Double) -> Double) -> Double {
let h = (to - from) / Double(n)
var x = from
var sum = 0.0
while x <= to - h {
sum += f(x + h)
x += h
}
return h * sum
}
private func integrateRectMid(from: Double, to: Double, n: Int, f: (Double) -> Double) -> Double {
let h = (to - from) / Double(n)
var x = from
var sum = 0.0
while x <= to - h {
sum += f(x + h / 2.0)
x += h
}
return h * sum
}
private func integrateTrapezium(from: Double, to: Double, n: Int, f: (Double) -> Double) -> Double {
let h = (to - from) / Double(n)
var sum = f(from) + f(to)
for i in 1..<n {
sum += 2 * f(from + Double(i) * h)
}
return h * sum / 2
}
private func integrateSimpson(from: Double, to: Double, n: Int, f: (Double) -> Double) -> Double {
let h = (to - from) / Double(n)
var sum1 = 0.0
var sum2 = 0.0
for i in 0..<n {
sum1 += f(from + h * Double(i) + h / 2.0)
}
for i in 1..<n {
sum2 += f(from + h * Double(i))
}
return h / 6.0 * (f(from) + f(to) + 4.0 * sum1 + 2.0 * sum2)
}
let types = IntegrationType.allCases
print("f(x) = x^3:", types.map({ integrate(from: 0, to: 1, n: 100, using: $0, f: { pow($0, 3) }) }))
print("f(x) = 1 / x:", types.map({ integrate(from: 1, to: 100, n: 1000, using: $0, f: { 1 / $0 }) }))
print("f(x) = x, 0 -> 5_000:", types.map({ integrate(from: 0, to: 5_000, n: 5_000_000, using: $0, f: { $0 }) }))
print("f(x) = x, 0 -> 6_000:", types.map({ integrate(from: 0, to: 6_000, n: 6_000_000, using: $0, f: { $0 }) }))</syntaxhighlight>
{{out}}
<pre>f(x) = x^3: [0.2450250000000004, 0.23532201000000041, 0.2401367512500004, 0.25002500000000005, 0.25000000000000006]
f(x) = 1 / x: [4.55599105751469, 4.654000076443428, 4.603772058385689, 4.60598605751468, 4.605170384957145]
f(x) = x, 0 -> 5_000: [12499997.500728704, 12499992.500729704, 12499995.000729209, 12500000.000000002, 12500000.0]
f(x) = x, 0 -> 6_000: [17999997.001390498, 17999991.0013915, 17999994.001391016, 18000000.000000004, 17999999.999999993]</pre>
=={{header|Tcl}}==
<
proc leftrect {f left right} {
Line 3,434 ⟶ 6,026:
puts [format " %-10s %s\t(%.1f%%)" $method $int $diff]
}
}</
<pre>integral of square(x) from 0 to 3.141592653589793 in 10 steps
leftrect 8.836788853885448 (-14.5%)
Line 3,448 ⟶ 6,040:
simpson 2.0000067844418012 (0.0%)</pre>
==
<!-- Not put into category so that it's considered unimplemented. -->
Line 3,464 ⟶ 6,056:
the integrand <math>f</math>, the bounds <math>(a,b)</math>, and the number of intervals <math>n</math>.
<
#import nat
#import flo
Line 3,470 ⟶ 6,062:
(integral_by "m") ("f","a","b","n") =
iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"</
An alternative way of defining this function shown below prevents redundant evaluations of the integrand
at the cost of building a table-driven finite map in advance.
<
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</
As mentioned in the Haskell solution, the latter choice is preferable if evaluating the integrand
is expensive.
An integrating function is defined for each method as follows.
<
right = integral_by "f". ("l","r"). "f" "r"
midpoint = integral_by "f". ("l","r"). "f" div\2. plus/"l" "r"
trapezium = integral_by "f". ("l","r"). div\2. plus "f"~~/"l" "r"
simpson = integral_by "f". ("l","r"). div\6. plus:-0. <"f" "l",times/4. "f" div\2. plus/"l" "r","f" "r"></
As shown above, the method passed to the <code>integral_by</code> function
is itself a higher order function taking an integrand <math>f</math> as an argument and
Line 3,489 ⟶ 6,081:
Here is a test program showing the results of integrating the square from zero to <math>\pi</math> in ten intervals
by all five methods.
<
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</
output:
<pre>
Line 3,503 ⟶ 6,095:
(The GNU Scientific Library integration routines are also callable in Ursala, and
are faster and more accurate.)
=={{header|VBA}}==
The following program does not follow the task requirement on two points: first, the same function is used for all quadrature methods, as they are really the same thing with different parameters (abscissas and weights). And since it's getting rather slow for a large number of intervals, the last two are integrated with resp. 50,000 and 60,000 intervals. It does not make sense anyway to use more, for such a simple function (and if really it were difficult to integrate, one would rely one more sophistcated methods).
<syntaxhighlight lang="vb">Option Explicit
Option Base 1
Function Quad(ByVal f As String, ByVal a As Double, _
ByVal b As Double, ByVal n As Long, _
ByVal u As Variant, ByVal v As Variant) As Double
Dim m As Long, h As Double, x As Double, s As Double, i As Long, j As Long
m = UBound(u)
h = (b - a) / n
s = 0#
For i = 1 To n
x = a + (i - 1) * h
For j = 1 To m
s = s + v(j) * Application.Run(f, x + h * u(j))
Next
Next
Quad = s * h
End Function
Function f1fun(x As Double) As Double
f1fun = x ^ 3
End Function
Function f2fun(x As Double) As Double
f2fun = 1 / x
End Function
Function f3fun(x As Double) As Double
f3fun = x
End Function
Sub Test()
Dim fun, f, coef, c
Dim i As Long, j As Long, s As Double
fun = Array(Array("f1fun", 0, 1, 100, 1 / 4), _
Array("f2fun", 1, 100, 1000, Log(100)), _
Array("f3fun", 0, 5000, 50000, 5000 ^ 2 / 2), _
Array("f3fun", 0, 6000, 60000, 6000 ^ 2 / 2))
coef = Array(Array("Left rect. ", Array(0, 1), Array(1, 0)), _
Array("Right rect. ", Array(0, 1), Array(0, 1)), _
Array("Midpoint ", Array(0.5), Array(1)), _
Array("Trapez. ", Array(0, 1), Array(0.5, 0.5)), _
Array("Simpson ", Array(0, 0.5, 1), Array(1 / 6, 4 / 6, 1 / 6)))
For i = 1 To UBound(fun)
f = fun(i)
Debug.Print f(1)
For j = 1 To UBound(coef)
c = coef(j)
s = Quad(f(1), f(2), f(3), f(4), c(2), c(3))
Debug.Print " " + c(1) + ": ", s, (s - f(5)) / f(5)
Next j
Next i
End Sub</syntaxhighlight>
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
var integrate = Fn.new { |a, b, n, f|
var h = (b - a) / n
var sum = List.filled(5, 0)
for (i in 0...n) {
var x = a + i * h
sum[0] = sum[0] + f.call(x)
sum[1] = sum[1] + f.call(x + h/2)
sum[2] = sum[2] + f.call(x + h)
sum[3] = sum[3] + (f.call(x) + f.call(x+h))/2
sum[4] = sum[4] + (f.call(x) + 4 * f.call(x + h/2) + f.call(x + h))/6
}
var methods = ["LeftRect ", "MidRect ", "RightRect", "Trapezium", "Simpson "]
for (i in 0..4) Fmt.print("$s = $h", methods[i], sum[i] * h)
System.print()
}
integrate.call(0, 1, 100) { |v| v * v * v }
integrate.call(1, 100, 1000) { |v| 1 / v }
integrate.call(0, 5000, 5000000) { |v| v }
integrate.call(0, 6000, 6000000) { |v| v }
</syntaxhighlight>
{{out}}
<pre>
LeftRect = 0.245025
MidRect = 0.249988
RightRect = 0.255025
Trapezium = 0.250025
Simpson = 0.25
LeftRect = 4.654991
MidRect = 4.604763
RightRect = 4.556981
Trapezium = 4.605986
Simpson = 4.60517
LeftRect = 12499997.5
MidRect = 12500000
RightRect = 12500002.5
Trapezium = 12500000
Simpson = 12500000
LeftRect = 17999997
MidRect = 18000000
RightRect = 18000003
Trapezium = 18000000
Simpson = 18000000
</pre>
=={{header|XPL0}}==
<syntaxhighlight lang="xpl0">include c:\cxpl\codes; \intrinsic 'code' declarations
func real Func(FN, X); \Return F(X) for function number FN
int FN; real X;
[case FN of
1: return X*X*X;
2: return 1.0/X;
3: return X
other return 0.0;
];
func Integrate(A, B, FN, N); \Display area under curve for function FN
real A, B; int FN, N; \limits A, B, and number of slices N
real DX, X, Area; \delta X
int I;
[DX:= (B-A)/float(N);
X:= A; Area:= 0.0; \rectangular left
for I:= 1 to N do
[Area:= Area + Func(FN,X)*DX; X:= X+DX];
RlOut(0, Area);
X:= A; Area:= 0.0; \rectangular right
for I:= 1 to N do
[X:= X+DX; Area:= Area + Func(FN,X)*DX];
RlOut(0, Area);
X:= A+DX/2.0; Area:= 0.0; \rectangular mid point
for I:= 1 to N do
[Area:= Area + Func(FN,X)*DX; X:= X+DX];
RlOut(0, Area);
X:= A; Area:= 0.0; \trapezium
for I:= 1 to N do
[Area:= Area + (Func(FN,X)+Func(FN,X+DX))/2.0*DX; X:= X+DX];
RlOut(0, Area);
X:= A; Area:= 0.0; \Simpson's rule
for I:= 1 to N do
[Area:= Area +
DX/6.0*(Func(FN,X) + 4.0*Func(FN,(X+X+DX)/2.0) + Func(FN,X+DX));
X:= X+DX];
RlOut(0, Area);
CrLf(0);
];
[Format(9,6);
Integrate(0.0, 1.0, 1, 100);
Integrate(1.0, 100.0, 2, 1000);
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
]</syntaxhighlight>
Interestingly, the small rounding errors creep in when millions of
approximations are done. If the five and six millions are changed to five
and six thousands then the rounding errors disappear. (They could have
been hidden by using scientific notation for the output format.)
{{out}}
<pre>
0.245025 0.255025 0.249988 0.250025 0.250000
4.654991 4.556981 4.604763 4.605986 4.605170
12499997.500729 12500002.500729 12500000.000729 12500000.000729 12500000.000729
17999997.001391 18000003.001391 18000000.001391 18000000.001391 18000000.001391
</pre>
=={{header|Yabasic}}==
Based on the XPL0entry and the Free BASIC entry
<syntaxhighlight lang="yabasic">// Rosetta Code problem: https://rosettacode.org/wiki/Numerical_integration
// by Jjuanhdez, 06/2022
print "function range steps leftrect midrect rightrect trap simpson "
frmt$ = "%1.10f"
print "f(x) = x^3 0 - 1 100 ";
Integrate(0.0, 1.0, 1, 100)
print "f(x) = 1/x 1 - 100 1000 ";
Integrate(1.0, 100.0, 2, 1000)
frmt$ = "%8.3f"
print "f(x) = x 0 - 5000 5000000 ";
Integrate(0.0, 5000.0, 3, 5000000)
print "f(x) = x 0 - 6000 6000000 ";
Integrate(0.0, 6000.0, 3, 6000000)
end
sub Func(FN, X) //Return F(X) for function number FN
switch FN
case 1
return X ^ 3
case 2
return 1.0 / X
case 3
return X
default
return 0.0
end switch
end sub
sub Integrate(A, B, FN, N) //Display area under curve for function FN
// A, B, FN limits A, B, and number of slices N
DX = (B-A)/N
X = A
Area = 0.0 //rectangular left
for i = 1 to N
Area = Area + Func(FN,X)*DX
X = X + DX
next i
print str$(Area, frmt$);
X = A
Area = 0.0 //rectangular right
for i = 1 to N
X = X + DX
Area = Area + Func(FN,X)*DX
next i
print " ";
print str$(Area, frmt$);
X = A + DX / 2.0
Area = 0.0 //rectangular mid point
for i = 1 to N
Area = Area + Func(FN,X)*DX
X = X + DX
next i
print " ";
print str$(Area, frmt$);
X = A
Area = 0.0 //trapezium
for i = 1 to N
Area = Area + (Func(FN,X)+Func(FN,X + DX))/2.0*DX
X = X + DX
next i
print " ";
print str$(Area, frmt$);
X = A
Area = 0.0 //Simpson's rule
for i = 1 to N
Area = Area + DX/6.0*(Func(FN,X) + 4.0*Func(FN,(X+X + DX)/2.0) + Func(FN,X + DX))
X = X + DX
next i
print " ";
print str$(Area, frmt$)
end sub</syntaxhighlight>
=={{header|zkl}}==
{{trans|D}}
<syntaxhighlight lang="zkl">fcn integrate(F,f,a,b,steps){
h:=(b - a) / steps;
h*(0).reduce(steps,'wrap(s,i){ F(f, h*i + a, h) + s },0.0);
}
fcn rectangularLeft(f,x) { f(x) }
fcn rectangularMiddle(f,x,h){ f(x+h/2) }
fcn rectangularRight(f,x,h) { f(x+h) }
fcn trapezium(f,x,h) { (f(x) + f(x+h))/2 }
fcn simpson(f,x,h) { (f(x) + 4.0*f(x+h/2) + f(x+h))/6 }
args:=T( T(fcn(x){ x.pow(3) }, 0.0, 1.0, 10),
T(fcn(x){ 1.0 / x }, 1.0, 100.0, 1000),
T(fcn(x){ x }, 0.0, 5000.0, 0d5_000_000),
T(fcn(x){ x }, 0.0, 6000.0, 0d6_000_000) );
fs:=T(rectangularLeft,rectangularMiddle,rectangularRight,
trapezium,simpson);
names:=fs.pump(List,"name",'+(":"),"%-18s".fmt);
foreach a in (args){
names.zipWith('wrap(nm,f){
"%s %f".fmt(nm,integrate(f,a.xplode())).println() }, fs);
println();
}</syntaxhighlight>
{{out}}
<pre>
rectangularLeft: 0.202500
rectangularMiddle: 0.248750
rectangularRight: 0.302500
trapezium: 0.252500
simpson: 0.250000
rectangularLeft: 4.654991
rectangularMiddle: 4.604763
rectangularRight: 4.556981
trapezium: 4.605986
simpson: 4.605170
rectangularLeft: 12499997.500000
rectangularMiddle: 12500000.000000
rectangularRight: 12500002.500000
trapezium: 12500000.000000
simpson: 12500000.000000
rectangularLeft: 17999997.000000
rectangularMiddle: 18000000.000000
rectangularRight: 18000003.000000
trapezium: 18000000.000000
simpson: 18000000.000000
</pre>
{{omit from|GUISS}}
|