Numerical integration: Difference between revisions
Content added Content deleted
(→{{header|PL/I}}: Rewritten to incorporate all the requirements.) |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 108: | Line 108: | ||
trace(simpson(f1, -1, 2 ,4 )); |
trace(simpson(f1, -1, 2 ,4 )); |
||
</lang> |
</lang> |
||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
Specification of a generic package implementing the five specified kinds of numerical integration: |
Specification of a generic package implementing the five specified kinds of numerical integration: |
||
Line 455: | Line 456: | ||
fun(x) { ; linear test function |
fun(x) { ; linear test function |
||
Return x |
Return x |
||
}</lang> |
}</lang> |
||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
||
Line 866: | Line 867: | ||
18000000 |
18000000 |
||
18000000</lang> |
18000000</lang> |
||
=={{header|C++}}== |
=={{header|C++}}== |
||
Line 1,095: | Line 1,097: | ||
simpsonsIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 0.0 |
simpsonsIntegration: calculated = 1.8e+07; exact = 1.8e+07; difference = 0.0 |
||
</lang> |
</lang> |
||
=={{header|CoffeeScript}}== |
=={{header|CoffeeScript}}== |
||
{{trans|python}} |
{{trans|python}} |
||
Line 1,834: | Line 1,837: | ||
end module FunctionHolder</lang> |
end module FunctionHolder</lang> |
||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
Based on the BASIC entry and the BBC BASIC entry |
Based on the BASIC entry and the BBC BASIC entry |
||
Line 2,977: | Line 2,981: | ||
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014 |
print integrate "i.trapezium "fn2 4 -1 2 ; 2.351014 |
||
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</lang> |
print integrate "i.simpsons "fn2 4 -1 2 ; 2.447732</lang> |
||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
<lang lua>function leftRect( f, a, b, n ) |
<lang lua>function leftRect( f, a, b, n ) |
||
Line 3,558: | Line 3,563: | ||
$_ |
$_ |
||
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|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}} |
|||
<lang perl6>use MONKEY-SEE-NO-EVAL; |
|||
sub leftrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a, $a+$h ... $b-$h; |
|||
} |
|||
sub rightrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a+$h, $a+$h+$h ... $b; |
|||
} |
|||
sub midrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a+$h/2, $a+$h+$h/2 ... $b-$h/2; |
|||
} |
|||
sub trapez(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
my $partial-sum += f($_) * 2 for $a+$h, $a+$h+$h ... $b-$h; |
|||
$h / 2 * [+] f($a), f($b), $partial-sum; |
|||
} |
|||
sub simpsons(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
my $h2 = $h/2; |
|||
my $sum1 = f($a + $h2); |
|||
my $sum2 = 0; |
|||
for $a+$h, *+$h ... $b-$h { |
|||
$sum1 += f($_ + $h2); |
|||
$sum2 += f($_); |
|||
} |
|||
($h / 6) * (f($a) + f($b) + 4*$sum1 + 2*$sum2); |
|||
} |
|||
sub integrate($f, $a, $b, $n, $exact) { |
|||
my @r0; |
|||
my $e = 0.000001; |
|||
@r0.push: "$f\n in [$a..$b] / $n\n"; |
|||
@r0.push: ' exact result: '~ $exact.round($e); |
|||
my (@r1,@r2,@r3,@r4,@r5); |
|||
my &f; |
|||
EVAL "&f = $f"; |
|||
my $p1 = Promise.start( { @r1.push: ' rectangle method left: '~ leftrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p2 = Promise.start( { @r2.push: ' rectangle method right: '~ rightrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p3 = Promise.start( { @r3.push: ' rectangle method mid: '~ midrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p4 = Promise.start( { @r4.push: 'composite trapezoidal rule: '~ trapez(&f, $a, $b, $n).round($e) } ); |
|||
my $p5 = Promise.start( { @r5.push: ' 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;</lang> |
|||
{{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 |
in [0..6000] / 6000000 |
||
exact result: 18000000 |
exact result: 18000000 |
||
Line 3,751: | Line 3,654: | ||
x 0 - 6000 6000000 17999997 18000000 18000003 18000000 18000000 |
x 0 - 6000 6000000 17999997 18000000 18000003 18000000 18000000 |
||
</pre> |
</pre> |
||
=={{header|PicoLisp}}== |
|||
<lang PicoLisp>(scl 6) |
|||
(de leftRect (Fun X) |
|||
(Fun X) ) |
|||
(de rightRect (Fun X H) |
|||
(Fun (+ X H)) ) |
|||
(de midRect (Fun X H) |
|||
(Fun (+ X (/ H 2))) ) |
|||
(de trapezium (Fun X H) |
|||
(/ (+ (Fun X) (Fun (+ X H))) 2) ) |
|||
(de simpson (Fun X H) |
|||
(*/ |
|||
(+ |
|||
(Fun X) |
|||
(* 4 (Fun (+ X (/ H 2)))) |
|||
(Fun (+ X H)) ) |
|||
6 ) ) |
|||
(de square (X) |
|||
(*/ X X 1.0) ) |
|||
(de integrate (Fun From To Steps Meth) |
|||
(let (H (/ (- To From) Steps) Sum 0) |
|||
(for (X From (>= (- To H) X) (+ X H)) |
|||
(inc 'Sum (Meth Fun X H)) ) |
|||
(*/ H Sum 1.0) ) ) |
|||
(prinl (round (integrate square 3.0 7.0 30 simpson)))</lang> |
|||
Output: |
|||
<pre>105.333</pre> |
|||
=={{header|PL/I}}== |
=={{header|PL/I}}== |
||
Line 3,832: | Line 3,771: | ||
1.799999700000000E+0007 1.800000000000000E+0007 1.800000300000000E+0007 1.800000000000000E+0007 1.800000000000000E+0007 |
1.799999700000000E+0007 1.800000000000000E+0007 1.800000300000000E+0007 1.800000000000000E+0007 1.800000000000000E+0007 |
||
</pre> |
</pre> |
||
=={{header|PicoLisp}}== |
|||
<lang PicoLisp>(scl 6) |
|||
(de leftRect (Fun X) |
|||
(Fun X) ) |
|||
(de rightRect (Fun X H) |
|||
(Fun (+ X H)) ) |
|||
(de midRect (Fun X H) |
|||
(Fun (+ X (/ H 2))) ) |
|||
(de trapezium (Fun X H) |
|||
(/ (+ (Fun X) (Fun (+ X H))) 2) ) |
|||
(de simpson (Fun X H) |
|||
(*/ |
|||
(+ |
|||
(Fun X) |
|||
(* 4 (Fun (+ X (/ H 2)))) |
|||
(Fun (+ X H)) ) |
|||
6 ) ) |
|||
(de square (X) |
|||
(*/ X X 1.0) ) |
|||
(de integrate (Fun From To Steps Meth) |
|||
(let (H (/ (- To From) Steps) Sum 0) |
|||
(for (X From (>= (- To H) X) (+ X H)) |
|||
(inc 'Sum (Meth Fun X H)) ) |
|||
(*/ H Sum 1.0) ) ) |
|||
(prinl (round (integrate square 3.0 7.0 30 simpson)))</lang> |
|||
Output: |
|||
<pre>105.333</pre> |
|||
=={{header|PureBasic}}== |
=={{header|PureBasic}}== |
||
Line 4,268: | Line 4,171: | ||
simpson: 17999999.999999993 |
simpson: 17999999.999999993 |
||
</lang> |
</lang> |
||
=={{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}} |
|||
<lang perl6>use MONKEY-SEE-NO-EVAL; |
|||
sub leftrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a, $a+$h ... $b-$h; |
|||
} |
|||
sub rightrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a+$h, $a+$h+$h ... $b; |
|||
} |
|||
sub midrect(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
$h * [+] do f($_) for $a+$h/2, $a+$h+$h/2 ... $b-$h/2; |
|||
} |
|||
sub trapez(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
my $partial-sum += f($_) * 2 for $a+$h, $a+$h+$h ... $b-$h; |
|||
$h / 2 * [+] f($a), f($b), $partial-sum; |
|||
} |
|||
sub simpsons(&f, $a, $b, $n) { |
|||
my $h = ($b - $a) / $n; |
|||
my $h2 = $h/2; |
|||
my $sum1 = f($a + $h2); |
|||
my $sum2 = 0; |
|||
for $a+$h, *+$h ... $b-$h { |
|||
$sum1 += f($_ + $h2); |
|||
$sum2 += f($_); |
|||
} |
|||
($h / 6) * (f($a) + f($b) + 4*$sum1 + 2*$sum2); |
|||
} |
|||
sub integrate($f, $a, $b, $n, $exact) { |
|||
my @r0; |
|||
my $e = 0.000001; |
|||
@r0.push: "$f\n in [$a..$b] / $n\n"; |
|||
@r0.push: ' exact result: '~ $exact.round($e); |
|||
my (@r1,@r2,@r3,@r4,@r5); |
|||
my &f; |
|||
EVAL "&f = $f"; |
|||
my $p1 = Promise.start( { @r1.push: ' rectangle method left: '~ leftrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p2 = Promise.start( { @r2.push: ' rectangle method right: '~ rightrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p3 = Promise.start( { @r3.push: ' rectangle method mid: '~ midrect(&f, $a, $b, $n).round($e) } ); |
|||
my $p4 = Promise.start( { @r4.push: 'composite trapezoidal rule: '~ trapez(&f, $a, $b, $n).round($e) } ); |
|||
my $p5 = Promise.start( { @r5.push: ' 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;</lang> |
|||
{{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}}== |
=={{header|REXX}}== |
||
Line 4,519: | Line 4,525: | ||
trapezium 1.98352353750945 (-0.8%) |
trapezium 1.98352353750945 (-0.8%) |
||
simpson 2.0000067844418 (0.0%)</pre> |
simpson 2.0000067844418 (0.0%)</pre> |
||
=={{header|Rust}}== |
=={{header|Rust}}== |
||
Line 4,634: | Line 4,639: | ||
(define t (integrate square 0 1 10 trapezium)) |
(define t (integrate square 0 1 10 trapezium)) |
||
(define s (integrate square 0 1 10 simpson))</lang> |
(define s (integrate square 0 1 10 simpson))</lang> |
||
=={{header|Sidef}}== |
|||
{{trans|Perl 6}} |
|||
<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);</lang> |
|||
=={{header|SequenceL}}== |
=={{header|SequenceL}}== |
||
Line 4,762: | Line 4,710: | ||
x 0 - 6000 17999997. 18000003. 18000000. 18000000. 18000000." |
x 0 - 6000 17999997. 18000003. 18000000. 18000000. 18000000." |
||
</pre> |
</pre> |
||
=={{header|Sidef}}== |
|||
{{trans|Perl 6}} |
|||
<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);</lang> |
|||
=={{header|Standard ML}}== |
=={{header|Standard ML}}== |