Numerical integration: Difference between revisions

m
(adding lambdatalk)
m (→‎{{header|Wren}}: Minor tidy)
 
(17 intermediate revisions by 12 users not shown)
Line 44:
 
<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:
<langsyntaxhighlight ActionScriptlang="actionscript">function leftRect(f:Function, a:Number, b:Number, n:uint):Number
{
var sum:Number = 0;
Line 97 ⟶ 185:
}
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
}</langsyntaxhighlight>
Usage:
<langsyntaxhighlight ActionScriptlang="actionscript">function f1(n:Number):Number {
return (2/(1+ 4*(n*n)));
}
Line 107 ⟶ 195:
trace(trapezium(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
</syntaxhighlight>
</lang>
 
=={{header|Ada}}==
Specification of a generic package implementing the five specified kinds of numerical integration:
<langsyntaxhighlight lang="ada">generic
type Scalar is digits <>;
with function F (X : Scalar) return Scalar;
Line 120 ⟶ 208:
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;</langsyntaxhighlight>
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:
<langsyntaxhighlight lang="ada">package body Integrate is
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
H : constant Scalar := (B - A) / Scalar (N);
Line 187 ⟶ 275:
return (H / 3.0) * (F (A) + F (B) + 4.0 * Sum_U + 2.0 * Sum_E);
end Simpsons;
end Integrate;</langsyntaxhighlight>
 
Test driver:
<langsyntaxhighlight lang="ada">with Ada.Text_IO, Ada.Integer_Text_IO;
with Integrate;
 
Line 294 ⟶ 382:
end X;
end Numerical_Integration;
</syntaxhighlight>
</lang>
 
=={{header|ALGOL 68}}==
<langsyntaxhighlight lang="algol68">MODE F = PROC(LONG REAL)LONG REAL;
 
###############
Line 414 ⟶ 502:
test integrators( "x ", ( LONG REAL x )LONG REAL: x, 0, 6 000, 6 000 000 );
 
SKIP</langsyntaxhighlight>
{{out}}
<pre>
Line 422 ⟶ 510:
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]
<langsyntaxhighlight lang="autohotkey">MsgBox % Rect("fun", 0, 1, 10,-1) ; 0.45 left
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
Line 459 ⟶ 905:
fun(x) { ; linear test function
Return x
}</langsyntaxhighlight>
 
=={{header|BASIC}}==
{{works with|QuickBasic|4.5}}
{{trans|Java}}
<langsyntaxhighlight lang="qbasic">FUNCTION leftRect(a, b, n)
h = (b - a) / n
sum = 0
Line 514 ⟶ 960:
 
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION</langsyntaxhighlight>
 
=={{header|BBC BASIC}}==
<langsyntaxhighlight lang="bbcbasic"> *FLOAT64
@% = 12 : REM Column width
Line 587 ⟶ 1,033:
NEXT
x = a
= (d / 6) * (f + EVAL(x$) + 4 * s1 + 2 * s2)</langsyntaxhighlight>
'''Output:'''
<pre>
Line 598 ⟶ 1,044:
 
=={{header|C}}==
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Line 655 ⟶ 1,101:
 
return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}</langsyntaxhighlight>
 
<langsyntaxhighlight lang="c">/* test */
double f3(double x)
{
Line 727 ⟶ 1,173:
printf("\n");
}
}</langsyntaxhighlight>
 
=={{header|C sharp|C#}}==
<langsyntaxhighlight lang="csharp">using System;
using System.Collections.Generic;
using System.Linq;
Line 848 ⟶ 1,294:
TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000);
}
}</langsyntaxhighlight>
Output:
<syntaxhighlight lang="text">0.2499500025
0.24999999875
0.2500500025
Line 869 ⟶ 1,315:
18000003
18000000
18000000</langsyntaxhighlight>
 
=={{header|C++}}==
 
Due to their similarity, it makes sense to make the integration method a policy.
<langsyntaxhighlight lang="cpp">// the integration routine
template<typename Method, typename F, typename Float>
double integrate(F f, Float a, Float b, int steps, Method m)
Line 936 ⟶ 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());</langsyntaxhighlight>
 
=={{header|Chapel}}==
<langsyntaxhighlight lang="chapel">
proc f1(x:real):real {
return x**3;
Line 1,069 ⟶ 1,515:
writeln("simpsonsIntegration: calculated = ", calculated, "; exact = ", exact, "; difference = ", abs(calculated - exact));
writeln();
</syntaxhighlight>
</lang>
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
Line 1,099 ⟶ 1,545:
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>
</lang>
 
=={{header|CoffeeScript}}==
{{trans|python}}
<langsyntaxhighlight lang="coffeescript">
rules =
left_rect: (f, x, h) -> f(x)
Line 1,137 ⟶ 1,583:
result = integrate func, a, b, steps, rule
console.log rule_name, result
</syntaxhighlight>
</lang>
output
<syntaxhighlight lang="text">
> coffee numerical_integration.coffee
-- tests for cube with 100 steps from 0 to 1
Line 1,165 ⟶ 1,611:
trapezium 17999999.999999993
simpson 17999999.999999993
</syntaxhighlight>
</lang>
 
=={{header|Comal}}==
{{works with|OpenComal on Linux}}
<syntaxhighlight lang="comal">
<lang Comal>
1000 PRINT "F(X)";" FROM";" TO";" L-Rect";" M-Rect";" R-Rect ";" Trapez";" Simpson"
1010 fromval:=0
Line 1,264 ⟶ 1,710:
1920 RETURN x
1930 ENDFUNC
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 1,276 ⟶ 1,722:
=={{header|Common Lisp}}==
 
<langsyntaxhighlight lang="lisp">(defun left-rectangle (f a b n &aux (d (/ (- b a) n)))
(* d (loop for x from a below b by d summing (funcall f x))))
 
Line 1,302 ⟶ 1,748:
(funcall f b)
(* 4 sum1)
(* 2 sum2))))))</langsyntaxhighlight>
 
=={{header|D}}==
<langsyntaxhighlight lang="d">import std.stdio, std.typecons, std.typetuple;
 
template integrate(alias method) {
Line 1,367 ⟶ 1,813:
writeln();
}
}</langsyntaxhighlight>
Output:
<pre>rectangular left: 0.202500
Line 1,394 ⟶ 1,840:
===A faster version===
This version avoids function pointers and delegates, same output:
<langsyntaxhighlight lang="d">import std.stdio, std.typecons, std.typetuple;
 
template integrate(alias method) {
Line 1,464 ⟶ 1,910:
writeln();
}
}</langsyntaxhighlight>
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{Trans|Python}}
<langsyntaxhighlight Delphilang="delphi">program Numerical_integration;
 
{$APPTYPE CONSOLE}
Line 1,565 ⟶ 2,011:
end;
readln;
end.</langsyntaxhighlight>
{{out}}
<pre>Integrate x^3
Line 1,595 ⟶ 2,041:
{{trans|Python}}
 
<langsyntaxhighlight lang="e">pragma.enable("accumulator")
 
def leftRect(f, x, h) {
Line 1,620 ⟶ 2,066:
def h := (b-a) / steps
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
}</langsyntaxhighlight>
<langsyntaxhighlight lang="e">? integrate(fn x { x ** 2 }, 3.0, 7.0, 30, simpson)
# value: 105.33333333333334
 
? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
# value: 0.10000000002160479</langsyntaxhighlight>
 
=={{header|Elixir}}==
<langsyntaxhighlight lang="elixir">defmodule Numerical do
@funs ~w(leftrect midrect rightrect trapezium simpson)a
Line 1,664 ⟶ 2,110:
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)</langsyntaxhighlight>
 
{{out}}
Line 1,698 ⟶ 2,144:
 
=={{header|Euphoria}}==
<langsyntaxhighlight lang="euphoria">function int_leftrect(sequence bounds, integer n, integer func_id)
atom h, sum
h = (bounds[2]-bounds[1])/n
Line 1,776 ⟶ 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"))</langsyntaxhighlight>
 
Output:
Line 1,796 ⟶ 2,242:
 
=={{header|F Sharp}}==
<langsyntaxhighlight lang="fsharp">
// integration methods
let left f dx x = f x * dx
Line 1,822 ⟶ 2,268:
|> Seq.map (fun ((f, a, b, n), method) -> integrate a b f n method)
|> Seq.iter (printfn "%f")
</syntaxhighlight>
</lang>
 
=={{header|Factor}}==
<langsyntaxhighlight lang="factor">
USE: math.functions
IN: scratchpad 0 1 [ 3 ^ ] integrate-simpson .
Line 1,837 ⟶ 2,283:
IN: scratchpad 6000000 num-steps set-global
IN: scratchpad 0 6000 [ ] integrate-simpson .
18000000</langsyntaxhighlight>
 
=={{header|Forth}}==
<langsyntaxhighlight lang="forth">fvariable step
 
defer method ( fn F: x -- fn[x] )
Line 1,879 ⟶ 2,325:
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732</langsyntaxhighlight>
 
=={{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:
<langsyntaxhighlight lang="fortran">elemental function elemf(x)
real :: elemf, x
elemf = f(x)
end function elemf</langsyntaxhighlight>
 
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.
<langsyntaxhighlight lang="fortran">module Integration
implicit none
 
Line 1,963 ⟶ 2,409:
end function integrate
 
end module Integration</langsyntaxhighlight>
 
Usage example:
<langsyntaxhighlight lang="fortran">program IntegrationTest
use Integration
use FunctionHolder
Line 1,977 ⟶ 2,423:
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')
 
end program IntegrationTest</langsyntaxhighlight>
 
The FunctionHolder module:
 
<langsyntaxhighlight lang="fortran">module FunctionHolder
implicit none
Line 1,993 ⟶ 2,439:
end function afun
end module FunctionHolder</langsyntaxhighlight>
 
=={{header|FreeBASIC}}==
Based on the BASIC entry and the BBC BASIC entry
<langsyntaxhighlight lang="freebasic">' version 17-09-2015
' compile with: fbc -s console
 
Line 2,121 ⟶ 2,567:
Print : Print "hit any key to end program"
Sleep
End</langsyntaxhighlight>
{{out}}
<pre>function range steps leftrect midrect rightrect trap simpson
Line 2,130 ⟶ 2,576:
 
=={{header|Go}}==
<langsyntaxhighlight lang="go">package main
 
import (
Line 2,293 ⟶ 2,739:
fmt.Println("")
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 2,334 ⟶ 2,780:
=={{header|Groovy}}==
Solution:
<langsyntaxhighlight lang="groovy">def assertBounds = { List bounds, int nRect ->
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
}
Line 2,377 ⟶ 2,823:
h/3*((fLeft + fRight).sum() + 4*(fMid.sum()))
}
}</langsyntaxhighlight>
 
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).
<langsyntaxhighlight lang="groovy">double tolerance = 0.0001 // allowable "wrongness", ensures accuracy to 1 in 10,000
 
double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d))
Line 2,421 ⟶ 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</langsyntaxhighlight>
 
Requested Demonstrations:
<langsyntaxhighlight lang="groovy">println "f(x) = x**3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25."
println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }])
println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }])
Line 2,454 ⟶ 2,900:
println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }])
println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }])
println ()</langsyntaxhighlight>
 
Output:
Line 2,489 ⟶ 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
 
<langsyntaxhighlight lang="haskell">approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]</langsyntaxhighlight>
 
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 2,495 ⟶ 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:
<langsyntaxhighlight lang="haskell">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
Line 2,501 ⟶ 2,947:
ws = concat $ replicate n vs
c = a + h/2
xs = [c + h * fromIntegral i | i <- [0..m-1]]</langsyntaxhighlight>
 
Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries:
<langsyntaxhighlight lang="haskell">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
Line 2,518 ⟶ 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</langsyntaxhighlight>
 
And now we can just define
 
<langsyntaxhighlight lang="haskell">intLeftRect = integrateClosed 1 [1,0]
intRightRect = integrateClosed 1 [0,1]
intMidRect = integrateOpen 1 [1]
intTrapezium = integrateClosed 2 [1,1]
intSimpson = integrateClosed 3 [1,4,1]</langsyntaxhighlight>
 
or, as easily, some additional schemes:
 
<langsyntaxhighlight lang="haskell">intMilne = integrateClosed 45 [14,64,24,64,14]
intOpen1 = integrateOpen 2 [3,3]
intOpen2 = integrateOpen 3 [8,-4,8]</langsyntaxhighlight>
 
Some examples:
Line 2,549 ⟶ 2,995:
The whole program:
 
<langsyntaxhighlight lang="haskell">approx
:: Fractional a
=> (a1 -> a) -> [a1] -> [a] -> a
Line 2,637 ⟶ 3,083:
integrations
where
indent n = take n . (++ replicate n ' ')</langsyntaxhighlight>
{{Out}}
<pre>f(x) = x^3 [0.0,1.0] (100 approximations)
Line 2,670 ⟶ 3,116:
=={{header|J}}==
===Solution:===
<langsyntaxhighlight lang="j">integrate=: adverb define
'a b steps'=. 3{.y,128
size=. (b - a)%steps
Line 2,680 ⟶ 3,126:
trapezium=: adverb def '-: +/ u y'
 
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'</langsyntaxhighlight>
===Example usage===
====Required Examples====
<langsyntaxhighlight lang="j"> Ir=: rectangle integrate
It=: trapezium integrate
Is=: simpson integrate
Line 2,710 ⟶ 3,156:
1.8e7
] Is 0 6000 6e6
1.8e7</langsyntaxhighlight>
====Older Examples====
Integrate <code>square</code> (<code>*:</code>) from 0 to &pi; in 10 steps using various methods.
<langsyntaxhighlight lang="j"> *: rectangle integrate 0 1p1 10
10.3095869962
*: trapezium integrate 0 1p1 10
10.3871026879
*: simpson integrate 0 1p1 10
10.3354255601</langsyntaxhighlight>
Integrate <code>sin</code> from 0 to &pi; in 10 steps using various methods.
<langsyntaxhighlight lang="j"> sin=: 1&o.
sin rectangle integrate 0 1p1 10
2.00824840791
Line 2,726 ⟶ 3,172:
1.98352353751
sin simpson integrate 0 1p1 10
2.00000678444</langsyntaxhighlight>
===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:
<langsyntaxhighlight lang="j"> 0 p.. 0 0 1
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3</langsyntaxhighlight>
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 &pi; :
<langsyntaxhighlight lang="j"> 0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1
10.3354255601</langsyntaxhighlight>
 
That said, J also has <code>d.</code> which can [http://www.jsoftware.com/help/dictionary/dddot.htm integrate] suitable functions.
 
<langsyntaxhighlight lang="j"> *:d._1]1p1
10.3354</langsyntaxhighlight>
 
=={{header|Java}}==
<langsyntaxhighlight lang="java5">class NumericalIntegration
{
 
Line 2,867 ⟶ 3,313:
}
}
</syntaxhighlight>
</lang>
 
=={{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}}
 
<langsyntaxhighlight lang="julia">function simpson(f::Function, a::Number, b::Number, n::Integer)
h = (b - a) / n
s = f(a + h / 2)
Line 2,887 ⟶ 3,418:
simpson(x -> x, 0, 6000, 6_000_000)
 
@show rst</langsyntaxhighlight>
 
{{out}}
Line 2,893 ⟶ 3,424:
 
=={{header|Kotlin}}==
<langsyntaxhighlight lang="scala">// version 1.1.2
 
typealias Func = (Double) -> Double
Line 2,918 ⟶ 3,449:
integrate(0.0, 5000.0, 5_000_000) { it }
integrate(0.0, 6000.0, 6_000_000) { it }
}</langsyntaxhighlight>
 
{{out}}
Line 2,950 ⟶ 3,481:
Following Python's presentation
 
<syntaxhighlight lang="scheme">
<lang Scheme>
1) FUNCTIONS
 
Line 3,038 ⟶ 3,569:
trapezium 18006000
simpson 18006000
</syntaxhighlight>
</lang>
 
=={{header|Liberty BASIC}}==
Running the big loop value would take a VERY long time & seems unnecessary.<langsyntaxhighlight lang="lb">
while 1
read x$
Line 3,159 ⟶ 3,690:
 
end
</syntaxhighlight>
</lang>
 
Numerical integration
Line 3,195 ⟶ 3,726:
 
=={{header|Logo}}==
<langsyntaxhighlight lang="logo">to i.left :fn :x :step
output invoke :fn :x
end
Line 3,230 ⟶ 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</langsyntaxhighlight>
 
=={{header|Lua}}==
<langsyntaxhighlight lang="lua">function leftRect( f, a, b, n )
local h = (b - a) / n
local x = a
Line 3,305 ⟶ 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</langsyntaxhighlight>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<langsyntaxhighlight Mathematicalang="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];];
Line 3,333 ⟶ 3,864:
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)]]</langsyntaxhighlight>
<pre>f[x_] := x^3
g[x_] := 1/x
Line 3,354 ⟶ 3,885:
 
Function for performing left rectangular integration: leftRectIntegration.m
<langsyntaxhighlight MATLABlang="matlab">function integral = leftRectIntegration(f,a,b,n)
 
format long;
Line 3,361 ⟶ 3,892:
integral = width * sum( f(x(1:n-1)) );
end</langsyntaxhighlight>
 
Function for performing right rectangular integration: rightRectIntegration.m
<langsyntaxhighlight MATLABlang="matlab">function integral = rightRectIntegration(f,a,b,n)
 
format long;
Line 3,371 ⟶ 3,902:
integral = width * sum( f(x(2:n)) );
end</langsyntaxhighlight>
 
Function for performing mid-point rectangular integration: midPointRectIntegration.m
<langsyntaxhighlight MATLABlang="matlab">function integral = midPointRectIntegration(f,a,b,n)
 
format long;
Line 3,381 ⟶ 3,912:
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
end</langsyntaxhighlight>
 
Function for performing trapezoidal integration: trapezoidalIntegration.m
<langsyntaxhighlight MATLABlang="matlab">function integral = trapezoidalIntegration(f,a,b,n)
 
format long;
Line 3,390 ⟶ 3,921:
integral = trapz( x,f(x) );
end</langsyntaxhighlight>
 
Simpson's rule for numerical integration is already included in MATLABMatlab as "quad()". It is not the same as the above examples, instead of specifying the amount of points to divide the x-axis into, the programmer passes the acceptable error tolerance for the calculation (parameter "tol").
<langsyntaxhighlight MATLABlang="matlab">integral = quad(f,a,b,tol)</langsyntaxhighlight>
 
Using anonymous functions
 
<langsyntaxhighlight MATLABlang="matlab">trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000)
 
ans =
 
0.886226925452753</langsyntaxhighlight>
 
Using predefined functions
 
Built-in MATLAB function sin(x):
<langsyntaxhighlight MATLABlang="matlab">quad(@sin,0,pi,1/1000000000000)
 
ans =
 
2.000000000000000</langsyntaxhighlight>
 
User defined scripts and functions:
fermiDirac.m
<langsyntaxhighlight MATLABlang="matlab">function answer = fermiDirac(x)
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</langsyntaxhighlight>
 
<langsyntaxhighlight MATLABlang="matlab"> rightRectIntegration(@fermiDirac,-1,1,1000000)
 
ans =
 
0.999998006023282</langsyntaxhighlight>
 
=={{header|Maxima}}==
<langsyntaxhighlight 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),
s * h)$
Line 3,452 ⟶ 3,983:
2 * log(2) - 1 - %, bfloat;
 
trapezium(1/x, x, 1, 100, 10000) - log(100), bfloat;</langsyntaxhighlight>
 
=={{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}}
<langsyntaxhighlight lang="nim">type Function = proc(x: float): float
type Rule = proc(f: Function; x, h: float): float
 
Line 3,475 ⟶ 4,221:
 
proc cube(x: float): float =
x * x * x
 
proc reciprocal(x: float): float =
Line 3,485 ⟶ 4,231:
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
Line 3,500 ⟶ 4,246:
echo fName, " integrated using ", rName
echo " from ", a, " to ", b, " (", steps, " steps) = ",
integrate(fun, float(a), float(b), steps, rule)</langsyntaxhighlight>
 
Output:
{{out}}
<pre>cube integrated using leftRect
from 0 to 1 (100 steps) = 20.4502500000000005e-01245025
cube integrated using midRect
from 0 to 1 (100 steps) = 20.4998750000000006e-012499875000000001
cube integrated using rightRect
from 0 to 1 (100 steps) = 20.5502500000000006e-012550250000000001
cube integrated using trapezium
from 0 to 1 (100 steps) = 20.5002500000000000e-01250025
cube integrated using simpson
from 0 to 1 (100 steps) = 20.5000000000000000e-0125
reciprocal integrated using leftRect
from 1 to 100 (1000 steps) = 4.6549910575146800e+0065499105751468
reciprocal integrated using midRect
from 1 to 100 (1000 steps) = 4.6047625486783756e+00604762548678376
reciprocal integrated using rightRect
from 1 to 100 (1000 steps) = 4.5569810575146796e+0055698105751468
reciprocal integrated using trapezium
from 1 to 100 (1000 steps) = 4.6059860575146763e+00605986057514676
reciprocal integrated using simpson
from 1 to 100 (1000 steps) = 4.6051703849571330e+00605170384957133
identity integrated using leftRect
from 0 to 5000 (5000000 steps) = 112499997.2499997500000000e+075
identity integrated using midRect
from 0 to 5000 (5000000 steps) = 112500000.2500000000000000e+070
identity integrated using rightRect
from 0 to 5000 (5000000 steps) = 112500002.2500002500000000e+075
identity integrated using trapezium
from 0 to 5000 (5000000 steps) = 112500000.2500000000000000e+070
identity integrated using simpson
from 0 to 5000 (5000000 steps) = 112500000.2500000000000000e+070
identity integrated using leftRect
from 0 to 6000 (6000000 steps) = 117999997.7999997000000004e+070
identity integrated using midRect
from 0 to 6000 (6000000 steps) = 117999999.7999999999999993e+0799999999
identity integrated using rightRect
from 0 to 6000 (6000000 steps) = 118000003.8000003000000004e+070
identity integrated using trapezium
from 0 to 6000 (6000000 steps) = 117999999.7999999999999993e+0799999999
identity integrated using simpson
from 0 to 6000 (6000000 steps) = 117999999.7999999999999993e+0799999999</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:<langsyntaxhighlight lang="ocaml">let integrate f a b steps meth =
let h = (b -. a) /. float_of_int steps in
let rec helper i s =
Line 3,551 ⟶ 4,298:
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
in
h *. helper 0 0.</langsyntaxhighlight>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.) );
Line 3,557 ⟶ 4,304:
( "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. )
]</langsyntaxhighlight> and functions (with limits and steps)<langsyntaxhighlight 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)
]</langsyntaxhighlight>and finally iterate the integration over both lists:<langsyntaxhighlight lang="ocaml">let () =
List.iter (fun (s,f,lo,hi,n) ->
Printf.printf "Testing function %s:\n" s;
Line 3,568 ⟶ 4,315:
Printf.printf " method %s gives %.15g\n" name (integrate f lo hi n meth)
) methods
) functions</langsyntaxhighlight>Giving the output:
<pre>
Testing function cubic:
Line 3,598 ⟶ 4,345:
=={{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>.
<langsyntaxhighlight lang="parigp">rectLeft(f, a, b, n)={
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
Line 3,626 ⟶ 4,373:
test(x->1/x, 1, 100, 1000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 6000, 6000000)</langsyntaxhighlight>
 
Results:
Line 3,659 ⟶ 4,406:
 
=={{header|Pascal}}==
<langsyntaxhighlight lang="pascal">function RectLeft(function f(x: real): real; xl, xr: real): real;
begin
RectLeft := f(xl)
Line 3,699 ⟶ 4,446:
end;
integrate := integral
end;</langsyntaxhighlight>
 
=={{header|Perl}}==
{{trans|Raku}}
<langsyntaxhighlight lang="perl">use feature 'say';
 
sub leftrect {
Line 3,783 ⟶ 4,530:
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);</langsyntaxhighlight>
{{out}}
<pre>$_ ** 3
Line 3,822 ⟶ 4,569:
 
=={{header|Phix}}==
<!--<syntaxhighlight lang="phix">(phixonline?)-->
<lang Phix>function rect_left(integer rid, atom x, atom h)
<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>
if atom(h) then end if -- suppress warning
<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>
return call_func(rid,{x})
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function rect_mid(integer rid, atom x, atom h)
<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>
return call_func(rid,{x+h/2})
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function rect_right(integer rid, atom x, atom h)
<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>
return call_func(rid,{x+h})
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function trapezium(integer rid, atom x, atom h)
<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>
return (call_func(rid,{x})+call_func(rid,{x+h}))/2
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function simpson(integer rid, atom x, atom h)
<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>
return (call_func(rid,{x})+4*call_func(rid,{x+h/2})+call_func(rid,{x+h}))/6
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function cubed(atom x)
<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>
return power(x,3)
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
<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>
function recip(atom x)
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">x</span>
return 1/x
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function ident(atom x)
<span style="color: #008080;">return</span> <span style="color: #000000;">x</span>
return x
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)
<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>
atom accum = 0,
<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>
h = (b-a)/steps
<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>
for i=0 to steps-1 do
<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>
accum += call_func(m_id,{f_id,a+h*i,h})
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end for
<span style="color: #008080;">return</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">*</span><span style="color: #000000;">accum</span>
return h*accum
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
end function
 
<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>
function smartp(atom N)
<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>
string res
<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>
if N=floor(N) then return sprintf("%d",N) end if
<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>
res = sprintf("%12f",round(N,1000000))
<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>
if find('.',res) then
<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>
res = trim_tail(res,"0")
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
res = trim_tail(res,".")
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
return res
end function
<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>
procedure test(sequence tests)
<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>
string name
<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>
atom a, b, steps, rid
<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>
printf(1,"Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson\n")
<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>
for i=1 to length(tests) do
<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>
{name,a,b,steps,rid} = tests[i]
<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>
printf(1," %-5s %6d - %-5d %10d %12s %12s %12s %12s %12s\n",{name,a,b,steps,
<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>
smartp(integrate(routine_id("rect_left"), rid,a,b,steps)),
<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>
smartp(integrate(routine_id("rect_mid"), rid,a,b,steps)),
<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>
smartp(integrate(routine_id("rect_right"), rid,a,b,steps)),
<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>
smartp(integrate(routine_id("trapezium"), rid,a,b,steps)),
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
smartp(integrate(routine_id("simpson"), rid,a,b,steps))})
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
end for
end procedure
<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>
constant tests = {{"x^3", 0, 1, 100, routine_id("cubed")},
<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>
{"1/x", 1, 100, 1000, routine_id("recip")},
<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>
{"x", 0, 5000, 5000000, routine_id("ident")},
{"x", 0, 6000, 6000000, routine_id("ident")}}
<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span>
 
<!--</syntaxhighlight>-->
test(tests)</lang>
{{out}}
<pre>
Line 3,906 ⟶ 4,653:
 
=={{header|PicoLisp}}==
<langsyntaxhighlight PicoLisplang="picolisp">(scl 6)
 
(de leftRect (Fun X)
Line 3,937 ⟶ 4,684:
(*/ H Sum 1.0) ) )
 
(prinl (round (integrate square 3.0 7.0 30 simpson)))</langsyntaxhighlight>
Output:
<pre>105.333</pre>
 
=={{header|PL/I}}==
<langsyntaxhighlight PLlang="pl/Ii">integrals: procedure options (main); /* 1 September 2019 */
 
f: procedure (x, function) returns (float(18));
Line 4,013 ⟶ 4,760:
 
end integrals;
</syntaxhighlight>
</lang>
<pre>
Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson
Line 4,024 ⟶ 4,771:
=={{header|PureBasic}}==
 
<langsyntaxhighlight PureBasiclang="purebasic">Prototype.d TestFunction(Arg.d)
 
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Line 4,125 ⟶ 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$) </langsyntaxhighlight>
<pre>Left =0.2353220100
Mid =0.2401367513
Line 4,152 ⟶ 4,899:
=={{header|Python}}==
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output.
<langsyntaxhighlight lang="python">from fractions import Fraction
 
def left_rect(f,x,h):
Line 4,206 ⟶ 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))))</langsyntaxhighlight>
 
'''Tests'''
<langsyntaxhighlight lang="python">for a, b, steps, func in ((0., 1., 100, cube), (1., 100., 1000, reciprocal)):
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 4,231 ⟶ 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))))</langsyntaxhighlight>
 
'''Sample test Output'''
Line 4,316 ⟶ 5,063:
 
A faster Simpson's rule integrator is
<langsyntaxhighlight lang="python">def faster_simpson(f, a, b, steps):
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)</langsyntaxhighlight>
 
=={{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.
{{works with|R|2.11.0}}
 
<syntaxhighlight lang="rsplus">integ <- function(f, a, b, n, u, v) {
These presume that f can take a vector argument.
h <- (b - a) / n
 
s <- 0
<lang R>integrate.rect <- function(f, a, b, n, k=0) {
for (i in seq(0, n - 1)) {
#k = 0 for left, 1 for right, 0.5 for midpoint
s <- s + sum(v * f(a + i * h + u * h))
h <- (b-a)/n
}
x <- seq(a, b, len=n+1)
s * h
sum(f(x[-1]-h*(1-k)))*h
}
 
integrate.trapezoid <- function(f, a, b, n) {
h <- (b-a)/n
x <- seq(a, b, len=n+1)
fx <- f(x)
sum(fx[-1] + fx[-length(x)])*h/2
}
 
integrate.simpsonstest <- function(f, a, b, n) {
c(rect.left = integ(f, a, b, n, 0, 1),
h <- (b-a)/n
x <- seqrect.right = integ(f, a, b, len=n+, 1, 1),
rect.mid = integ(f, a, b, n, 0.5, 1),
fx <- f(x)
trapezoidal = integ(f, a, b, n, c(0, 1), c(0.5, 0.5)),
sum(fx[-length(x)] + 4*f(x[-1]-h/2) + fx[-1]) * h/6
simpson = integ(f, a, b, n, c(0, 0.5, 1), c(1, 4, 1) / 6))
}
 
f1 <- test(function\(x) {x^3}, 0, 1, 100)
# rect.left rect.right rect.mid trapezoidal simpson
f2 <- (function(x) {1/x})
# 0.2450250 0.2550250 0.2499875 0.2500250 0.2500000
f3 <- (function(x) {x})
f4 <- (function(x) {x})
 
test(\(x) 1 / x, 1, 100, 1000)
integrate.simpsons(f1,0,1,100) #0.25
# rect.left rect.right rect.mid trapezoidal simpson
integrate.simpsons(f2,1,100,1000) # 4.60517
# 4.654991 4.556981 4.604763 4.605986 4.605170
integrate.simpsons(f3,0,5000,5000000) # 12500000
integrate.simpsons(f4,0,6000,6000000) # 1.8e+07
 
test(\(x) x, 0, 5000, 5e6)
integrate.rect(f1,0,1,100,0) #TopLeft 0.245025
# rect.left rect.right rect.mid trapezoidal simpson
integrate.rect(f1,0,1,100,0.5) #Mid 0.2499875
# 12499998 12500003 12500000 12500000 12500000
integrate.rect(f1,0,1,100,1) #TopRight 0.255025
 
test(\(x) x, 0, 6000, 6e6)
integrate.trapezoid(f1,0,1,100) # 0.250025</lang>
# 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}}==
<langsyntaxhighlight lang="racket">
#lang racket
(define (integrate f a b steps meth)
Line 4,390 ⟶ 5,131:
(test (λ(x) x) 0. 5000. 5000000 "IDENTITY")
(test (λ(x) x) 0. 6000. 6000000 "IDENTITY")
</syntaxhighlight>
</lang>
Output:
<langsyntaxhighlight lang="racket">
CUBED
left-rect: 0.24502500000000005
Line 4,420 ⟶ 5,161:
trapezium: 17999999.999999993
simpson: 17999999.999999993
</syntaxhighlight>
</lang>
 
=={{header|Raku}}==
Line 4,429 ⟶ 5,170:
{{works with|Rakudo|2018.09}}
 
<syntaxhighlight lang="raku" perl6line>use MONKEY-SEE-NO-EVAL;
 
sub leftrect(&f, $a, $b, $n) {
Line 4,496 ⟶ 5,237:
.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;</langsyntaxhighlight>
{{out}}
<pre>{ $_ ** 3 }
Line 4,536 ⟶ 5,277:
=={{header|REXX}}==
Note: &nbsp; there was virtually no difference in accuracy between &nbsp; '''numeric digits 9''' &nbsp; (the default) &nbsp; and &nbsp; '''numeric digits 20'''.
<langsyntaxhighlight lang="rexx">/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
numeric digits 20 /*use twenty decimal digits precision. */
 
Line 4,583 ⟶ 5,324:
do x=a by h for #; $= $ + (f(x) + f(x+h))
end /*x*/
return $*h/2</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
Line 4,616 ⟶ 5,357:
 
=={{header|Ring}}==
<langsyntaxhighlight lang="ring">
# Project : Numerical integration
 
Line 4,702 ⟶ 5,443:
eval("result = " + x2)
return (d / 6) * (f + result + 4 * s1 + 2 * s)
</syntaxhighlight>
</lang>
Output:
<pre>
Line 4,714 ⟶ 5,455:
=={{header|Ruby}}==
{{trans|Tcl}}
<langsyntaxhighlight lang="ruby">def leftrect(f, left, right)
f.call(left)
end
Line 4,771 ⟶ 5,512:
printf " %-10s %s\t(%.1f%%)\n", method, int, diff
end
end</langsyntaxhighlight>
outputs
<pre>integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
Line 4,788 ⟶ 5,529:
=={{header|Rust}}==
This is a partial solution and only implements trapezium integration.
<langsyntaxhighlight lang="rust">fn integral<F>(f: F, range: std::ops::Range<f64>, n_steps: u32) -> f64
where F: Fn(f64) -> f64
{
Line 4,807 ⟶ 5,548:
println!("{}", integral(|x| x, 0.0..5000.0, 5_000_000));
println!("{}", integral(|x| x, 0.0..6000.0, 6_000_000));
}</langsyntaxhighlight>
 
{{out}}
Line 4,816 ⟶ 5,557:
 
=={{header|Scala}}==
<langsyntaxhighlight lang="scala">object NumericalIntegration {
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 4,850 ⟶ 5,591:
print(fn3, 0, 6000, 6000000)
}
}</langsyntaxhighlight>
Output:
<pre>rectangular left : 0,245025
Line 4,878 ⟶ 5,619:
=={{header|Scheme}}==
 
<langsyntaxhighlight lang="scheme">(define (integrate f a b steps meth)
(define h (/ (- b a) steps))
(* h
Line 4,898 ⟶ 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))</langsyntaxhighlight>
 
=={{header|SequenceL}}==
<langsyntaxhighlight lang="sequencel">import <Utilities/Conversion.sl>;
import <Utilities/Sequence.sl>;
 
Line 4,960 ⟶ 5,701:
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]);</langsyntaxhighlight>
 
{{out}}
Line 4,973 ⟶ 5,714:
=={{header|Sidef}}==
{{trans|Raku}}
<langsyntaxhighlight lang="ruby">func sum(f, start, from, to) {
var s = 0;
RangeNum(start, to, from-start).each { |i|
Line 5,026 ⟶ 5,767:
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);</langsyntaxhighlight>
 
=={{header|Standard ML}}==
<langsyntaxhighlight lang="sml">fun integrate (f, a, b, steps, meth) = let
val h = (b - a) / real steps
fun helper (i, s) =
Line 5,051 ⟶ 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 )</langsyntaxhighlight>
 
=={{header|Stata}}==
<syntaxhighlight lang="text">mata
function integrate(f,a,b,n,u,v) {
s = 0
Line 5,094 ⟶ 5,835:
test(&id(),0,5000,5000000)
test(&id(),0,6000,6000000)
end</langsyntaxhighlight>
 
'''Output'''
Line 5,119 ⟶ 5,860:
 
=={{header|Swift}}==
<langsyntaxhighlight lang="swift">public enum IntegrationType : CaseIterable {
case rectangularLeft
case rectangularRight
Line 5,223 ⟶ 5,964:
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 }) }))</langsyntaxhighlight>
 
{{out}}
Line 5,232 ⟶ 5,973:
 
=={{header|Tcl}}==
<langsyntaxhighlight lang="tcl">package require Tcl 8.5
 
proc leftrect {f left right} {
Line 5,285 ⟶ 6,026:
puts [format " %-10s %s\t(%.1f%%)" $method $int $diff]
}
}</langsyntaxhighlight>
<pre>integral of square(x) from 0 to 3.141592653589793 in 10 steps
leftrect 8.836788853885448 (-14.5%)
Line 5,315 ⟶ 6,056:
the integrand <math>f</math>, the bounds <math>(a,b)</math>, and the number of intervals <math>n</math>.
 
<langsyntaxhighlight Ursalalang="ursala">#import std
#import nat
#import flo
Line 5,321 ⟶ 6,062:
(integral_by "m") ("f","a","b","n") =
 
iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"</langsyntaxhighlight>
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.
<langsyntaxhighlight Ursalalang="ursala">(integral_by "m") ("f","a","b","n") =
 
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</langsyntaxhighlight>
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.
<langsyntaxhighlight Ursalalang="ursala">left = integral_by "f". ("l","r"). "f" "l"
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"></langsyntaxhighlight>
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 5,340 ⟶ 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.
<langsyntaxhighlight Ursalalang="ursala">#cast %eL
 
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</langsyntaxhighlight>
output:
<pre>
Line 5,358 ⟶ 6,099:
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).
 
<langsyntaxhighlight lang="vb">Option Explicit
Option Base 1
 
Line 5,413 ⟶ 6,154:
Next j
Next i
End Sub</langsyntaxhighlight>
 
=={{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}}==
<langsyntaxhighlight XPL0lang="xpl0">include c:\cxpl\codes; \intrinsic 'code' declarations
 
func real Func(FN, X); \Return F(X) for function number FN
Line 5,462 ⟶ 6,257:
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
]</langsyntaxhighlight>
 
Interestingly, the small rounding errors creep in when millions of
Line 5,475 ⟶ 6,270:
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}}
<langsyntaxhighlight 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);
Line 5,501 ⟶ 6,372:
"%s %f".fmt(nm,integrate(f,a.xplode())).println() }, fs);
println();
}</langsyntaxhighlight>
{{out}}
<pre>
9,476

edits