Numerical integration/Adaptive Simpson's method: Difference between revisions
Content added Content deleted
m (→{{header|Lambdatalk}}: minor edit) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 36: | Line 36: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="11l">T Triple |
||
Float m, fm, simp |
Float m, fm, simp |
||
F (m, fm, simp) |
F (m, fm, simp) |
||
Line 75: | Line 75: | ||
V (a, b) = (0.0, 1.0) |
V (a, b) = (0.0, 1.0) |
||
V sinx = quad_asr(x -> sin(x), a, b, 1e-09) |
V sinx = quad_asr(x -> sin(x), a, b, 1e-09) |
||
print(‘Simpson's integration of sine from #. to #. = #.’.format(a, b, sinx))</ |
print(‘Simpson's integration of sine from #. to #. = #.’.format(a, b, sinx))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 84: | Line 84: | ||
=={{header|C}}== |
=={{header|C}}== |
||
{{trans|zkl}} |
{{trans|zkl}} |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 123: | Line 123: | ||
printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx); |
printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 132: | Line 132: | ||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{trans|Julia}} |
{{trans|Julia}} |
||
< |
<syntaxhighlight lang="factor">USING: formatting kernel locals math math.functions math.ranges |
||
sequences ; |
sequences ; |
||
IN: rosetta-code.simpsons |
IN: rosetta-code.simpsons |
||
Line 145: | Line 145: | ||
[ sin ] 0 1 100 simps |
[ sin ] 0 1 100 simps |
||
"Simpson's rule integration of sin from 0 to 1 is: %u\n" printf</ |
"Simpson's rule integration of sin from 0 to 1 is: %u\n" printf</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 153: | Line 153: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article. |
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 195: | Line 195: | ||
sinx := quadAsr(math.Sin, a, b, 1e-09) |
sinx := quadAsr(math.Sin, a, b, 1e-09) |
||
fmt.Printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx) |
fmt.Printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 218: | Line 218: | ||
<syntaxhighlight lang="j"> |
|||
<lang J> |
|||
Note'expected answer computed by j www.jsoftware.com' |
Note'expected answer computed by j www.jsoftware.com' |
||
Line 264: | Line 264: | ||
u uquad_asr a, fa, b, fb, eps, t Simp, t M, t Fm |
u uquad_asr a, fa, b, fb, eps, t Simp, t M, t Fm |
||
) |
) |
||
</syntaxhighlight> |
|||
</lang> |
|||
<pre> |
<pre> |
||
Line 272: | Line 272: | ||
=={{header|Java}}== |
=={{header|Java}}== |
||
<syntaxhighlight lang="java"> |
|||
<lang Java> |
|||
import java.util.function.Function; |
import java.util.function.Function; |
||
Line 325: | Line 325: | ||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 334: | Line 334: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Originally from Modesto Mas, https://mmas.github.io/simpson-integration-julia |
Originally from Modesto Mas, https://mmas.github.io/simpson-integration-julia |
||
< |
<syntaxhighlight lang="julia">function simps(f::Function, a::Number, b::Number, n::Number) |
||
iseven(n) || throw("n must be even, and $n was given") |
iseven(n) || throw("n must be even, and $n was given") |
||
h = (b-a)/n |
h = (b-a)/n |
||
Line 344: | Line 344: | ||
println("Simpson's rule integration of sin from 0 to 1 is: ", simps(sin, 0.0, 1.0, 100)) |
println("Simpson's rule integration of sin from 0 to 1 is: ", simps(sin, 0.0, 1.0, 100)) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936 |
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936 |
||
Line 351: | Line 351: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="scala">// Version 1.2.71 |
||
import kotlin.math.abs |
import kotlin.math.abs |
||
Line 393: | Line 393: | ||
val sinx = quadAsr(::sin, a, b, 1.0e-09) |
val sinx = quadAsr(::sin, a, b, 1.0e-09) |
||
println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}") |
println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}") |
||
}</ |
}</syntaxhighlight> |
||
{{output}} |
{{output}} |
||
Line 401: | Line 401: | ||
=={{header|Lambdatalk}}== |
=={{header|Lambdatalk}}== |
||
<syntaxhighlight lang="scheme"> |
|||
<lang Scheme> |
|||
1) using the sin's primitive -cos(x) |
1) using the sin's primitive -cos(x) |
||
Line 452: | Line 452: | ||
-> 0.45969769413186023 |
-> 0.45969769413186023 |
||
(0.45969769413186023) // reference value |
(0.45969769413186023) // reference value |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
Direct translation of Python code from Wikipedia entry. |
Direct translation of Python code from Wikipedia entry. |
||
< |
<syntaxhighlight lang="nim">import math, sugar |
||
type Func = float -> float |
type Func = float -> float |
||
Line 485: | Line 485: | ||
result = f.quadAsr(a, fa, b, fb, eps, whole, m, fm) |
result = f.quadAsr(a, fa, b, fb, eps, whole, m, fm) |
||
echo "Simpson's integration of sine from 0 to 1 = ", sin.quadAsr(0, 1, 1e-9)</ |
echo "Simpson's integration of sine from 0 to 1 = ", sin.quadAsr(0, 1, 1e-9)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 492: | Line 492: | ||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
Line 523: | Line 523: | ||
my ($a, $b) = (0, 1); |
my ($a, $b) = (0, 1); |
||
my $sin = adaptive_Simpson_quadrature('sin', $a, $b, 1e-9); |
my $sin = adaptive_Simpson_quadrature('sin', $a, $b, 1e-9); |
||
printf "Simpson's integration of sine from $a to $b = %.9f", $sin</ |
printf "Simpson's integration of sine from $a to $b = %.9f", $sin</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre> |
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre> |
||
Line 529: | Line 529: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #008080;">function</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</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: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">function</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</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: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">)</span> |
||
Line 569: | Line 569: | ||
<span style="color: #000000;">sinx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadAsr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">_sin</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;">1e-09</span><span style="color: #0000FF;">)</span> |
<span style="color: #000000;">sinx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadAsr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">_sin</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;">1e-09</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Simpson's integration of sine from %g to %g = %f\n"</span><span style="color: #0000FF;">,</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;">sinx</span><span style="color: #0000FF;">})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Simpson's integration of sine from %g to %g = %f\n"</span><span style="color: #0000FF;">,</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;">sinx</span><span style="color: #0000FF;">})</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 576: | Line 576: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
<syntaxhighlight lang="python"> |
|||
<lang Python> |
|||
#! python3 |
#! python3 |
||
Line 637: | Line 637: | ||
main() |
main() |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
Line 644: | Line 644: | ||
Fairly direct translation of the Python code. |
Fairly direct translation of the Python code. |
||
<lang |
<syntaxhighlight lang="raku" line>sub adaptive-Simpson-quadrature(&f, $left, $right, \ε = 1e-9) { |
||
my $lf = f($left); |
my $lf = f($left); |
||
my $rf = f($right); |
my $rf = f($right); |
||
Line 669: | Line 669: | ||
my ($a, $b) = 0e0, 1e0; |
my ($a, $b) = 0e0, 1e0; |
||
my $sin = adaptive-Simpson-quadrature(&sin, $a, $b, 1e-9).round(10**-9);; |
my $sin = adaptive-Simpson-quadrature(&sin, $a, $b, 1e-9).round(10**-9);; |
||
say "Simpson's integration of sine from $a to $b = $sin";</ |
say "Simpson's integration of sine from $a to $b = $sin";</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre> |
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre> |
||
Line 675: | Line 675: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="rexx">/*REXX program performs numerical integration using adaptive Simpson's method. */ |
||
numeric digits length( pi() ) - length(.) /*use # of digits in pi for precision. */ |
numeric digits length( pi() ) - length(.) /*use # of digits in pi for precision. */ |
||
a= 0; b= 1; f= 'SIN' /*define values for A, B, and F. */ |
a= 0; b= 1; f= 'SIN' /*define values for A, B, and F. */ |
||
Line 707: | Line 707: | ||
do k=2 by 2 until p=#; p= #; _= - _ * q / (k * (k+1) ); #= # + _ |
do k=2 by 2 until p=#; p= #; _= - _ * q / (k * (k+1) ); #= # + _ |
||
end /*k*/ |
end /*k*/ |
||
return #</ |
return #</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 715: | Line 715: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="ruby">func adaptive_Simpson_quadrature(f, left, right, ε = 1e-9) { |
||
func quadrature_mid(l, lf, r, rf) { |
func quadrature_mid(l, lf, r, rf) { |
||
Line 740: | Line 740: | ||
var (a = 0, b = 1) |
var (a = 0, b = 1) |
||
var area = adaptive_Simpson_quadrature({ .sin }, a, b, 1e-15).round(-15) |
var area = adaptive_Simpson_quadrature({ .sin }, a, b, 1e-15).round(-15) |
||
say "Simpson's integration of sine from #{a} to #{b} ≈ #{area}"</ |
say "Simpson's integration of sine from #{a} to #{b} ≈ #{area}"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 749: | Line 749: | ||
{{trans|Go}} |
{{trans|Go}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="ecmascript">import "/fmt" for Fmt |
||
/* "structured" adaptive version, translated from Racket */ |
/* "structured" adaptive version, translated from Racket */ |
||
Line 792: | Line 792: | ||
var b = 1 |
var b = 1 |
||
var sinx = quadAsr.call(Fn.new { |x| x.sin }, a, b, 1e-09) |
var sinx = quadAsr.call(Fn.new { |x| x.sin }, a, b, 1e-09) |
||
Fmt.print("Simpson's integration of sine from $d to $d = $f", a, b, sinx)</ |
Fmt.print("Simpson's integration of sine from $d to $d = $f", a, b, sinx)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 801: | Line 801: | ||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="zkl"># "structured" adaptive version, translated from Racket |
||
fcn _quad_simpsons_mem(f, a,fa, b,fb){ |
fcn _quad_simpsons_mem(f, a,fa, b,fb){ |
||
#Evaluates the Simpson's Rule, also returning m and f(m) to reuse""" |
#Evaluates the Simpson's Rule, also returning m and f(m) to reuse""" |
||
Line 823: | Line 823: | ||
m,fm,whole := _quad_simpsons_mem(f, a,fa, b,fb); |
m,fm,whole := _quad_simpsons_mem(f, a,fa, b,fb); |
||
_quad_asr(f, a,fa, b,fb, eps,whole,m,fm); |
_quad_asr(f, a,fa, b,fb, eps,whole,m,fm); |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="zkl">sinx:=quad_asr((1.0).sin.unbind(), 0.0, 1.0, 1e-09); |
||
println("Simpson's integration of sine from 1 to 2 = ",sinx);</ |
println("Simpson's integration of sine from 1 to 2 = ",sinx);</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Revision as of 00:13, 28 August 2022
Numerical integration/Adaptive Simpson's method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Lychee (1969)'s Modified Adaptive Simpson's method (doi:10.1145/321526.321537) is a numerical quadrature method that recursively bisects the interval until the precision is high enough.
; Lychee's ASR, Modifications 1, 2, 3 procedure _quad_asr_simpsons(f, a, fa, b, fb) m := (a + b) / 2 fm := f(m) h := b - a return multiple [m, fm, (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)] procedure _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, depth) lm, flm, left := _quad_asr_simpsons(f, a, fa, m, fm) rm, frm, right := _quad_asr_simpsons(f, m, fm, b, fb) delta := left + right - whole tol' := tol / 2 if depth <= 0 or tol' == tol or abs(delta) <= 15 * tol: return left + right + delta / 15 else: return _quad_asr(f, a, fa, m, fm, tol', left , lm, flm, depth - 1) + _quad_asr(f, m, fm, b, fb, tol', right, rm, frm, depth - 1) procedure quad_asr(f, a, b, tol, depth) fa := f(a) fb := f(b) m, fm, whole := _quad_asr_simpsons(f, a, fa, b, fb) return _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, depth) |
11l
T Triple
Float m, fm, simp
F (m, fm, simp)
.m = m
.fm = fm
.simp = simp
F _quad_simpsons_mem(f, Float a, fa, b, fb)
‘
Evaluates Simpson's Rule, also returning m and f(m) to reuse.
’
V m = a + (b - a) / 2
V fm = f(m)
V simp = abs(b - a) / 6 * (fa + 4 * fm + fb)
R Triple(m, fm, simp)
F _quad_asr(f, Float a, fa, b, fb, eps, whole, m, fm) -> Float
‘
Efficient recursive implementation of adaptive Simpson's rule.
Function values at the start, middle, end of the intervals are retained.
’
V lt = _quad_simpsons_mem(f, a, fa, m, fm)
V rt = _quad_simpsons_mem(f, m, fm, b, fb)
V delta = lt.simp + rt.simp - whole
R (I (abs(delta) <= eps * 15) {lt.simp + rt.simp + delta / 15}
E _quad_asr(f, a, fa, m, fm, eps / 2, lt.simp, lt.m, lt.fm)
+ _quad_asr(f, m, fm, b, fb, eps / 2, rt.simp, rt.m, rt.fm))
F quad_asr(f, Float a, b, eps) -> Float
‘
Integrate f from a to b using ASR with max error of eps.
’
V fa = f(a)
V fb = f(b)
V t = _quad_simpsons_mem(f, a, fa, b, fb)
R _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm)
V (a, b) = (0.0, 1.0)
V sinx = quad_asr(x -> sin(x), a, b, 1e-09)
print(‘Simpson's integration of sine from #. to #. = #.’.format(a, b, sinx))
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
C
#include <stdio.h>
#include <math.h>
typedef struct { double m; double fm; double simp; } triple;
/* "structured" adaptive version, translated from Racket */
triple _quad_simpsons_mem(double (*f)(double), double a, double fa, double b, double fb) {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse.
double m = (a + b) / 2;
double fm = f(m);
double simp = fabs(b - a) / 6 * (fa + 4*fm + fb);
triple t = {m, fm, simp};
return t;
}
double _quad_asr(double (*f)(double), double a, double fa, double b, double fb, double eps, double whole, double m, double fm) {
// Efficient recursive implementation of adaptive Simpson's rule.
// Function values at the start, middle, end of the intervals are retained.
triple lt = _quad_simpsons_mem(f, a, fa, m, fm);
triple rt = _quad_simpsons_mem(f, m, fm, b, fb);
double delta = lt.simp + rt.simp - whole;
if (fabs(delta) <= eps * 15) return lt.simp + rt.simp + delta/15;
return _quad_asr(f, a, fa, m, fm, eps/2, lt.simp, lt.m, lt.fm) +
_quad_asr(f, m, fm, b, fb, eps/2, rt.simp, rt.m, rt.fm);
}
double quad_asr(double (*f)(double), double a, double b, double eps) {
// Integrate f from a to b using ASR with max error of eps.
double fa = f(a);
double fb = f(b);
triple t = _quad_simpsons_mem(f, a, fa, b, fb);
return _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm);
}
int main(){
double a = 0.0, b = 1.0;
double sinx = quad_asr(sin, a, b, 1e-09);
printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx);
return 0;
}
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Factor
USING: formatting kernel locals math math.functions math.ranges
sequences ;
IN: rosetta-code.simpsons
:: simps ( f a b n -- x )
n even?
[ n "n must be even; %d was given" sprintf throw ] unless
b a - n / :> h
1 n 2 <range> 2 n 1 - 2 <range>
[ [ a + h * f call ] map-sum ] bi@ [ 4 ] [ 2 ] bi*
[ * ] 2bi@ a b [ f call ] bi@ + + + h 3 / * ; inline
[ sin ] 0 1 100 simps
"Simpson's rule integration of sin from 0 to 1 is: %u\n" printf
- Output:
Simpson's rule integration of sin from 0 to 1 is: 0.4596976941573994
Go
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article.
package main
import (
"fmt"
"math"
)
type F = func(float64) float64
/* "structured" adaptive version, translated from Racket */
func quadSimpsonsMem(f F, a, fa, b, fb float64) (m, fm, simp float64) {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse.
m = (a + b) / 2
fm = f(m)
simp = math.Abs(b-a) / 6 * (fa + 4*fm + fb)
return
}
func quadAsrRec(f F, a, fa, b, fb, eps, whole, m, fm float64) float64 {
// Efficient recursive implementation of adaptive Simpson's rule.
// Function values at the start, middle, end of the intervals are retained.
lm, flm, left := quadSimpsonsMem(f, a, fa, m, fm)
rm, frm, right := quadSimpsonsMem(f, m, fm, b, fb)
delta := left + right - whole
if math.Abs(delta) <= eps*15 {
return left + right + delta/15
}
return quadAsrRec(f, a, fa, m, fm, eps/2, left, lm, flm) +
quadAsrRec(f, m, fm, b, fb, eps/2, right, rm, frm)
}
func quadAsr(f F, a, b, eps float64) float64 {
// Integrate f from a to b using ASR with max error of eps.
fa, fb := f(a), f(b)
m, fm, whole := quadSimpsonsMem(f, a, fa, b, fb)
return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}
func main() {
a, b := 0.0, 1.0
sinx := quadAsr(math.Sin, a, b, 1e-09)
fmt.Printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx)
}
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
J
Typically one would choose the library implementation:
load'~addons/math/misc/integrat.ijs' NB. integrate returns definite integral and estimated digits of accuracy 1&o. integrate 0 1 0.459698 9 NB. adapt implements adaptive Simpson's rule, however recomputes some integrands 1&o. adapt 0 1 1e_9 0.459698
Note'expected answer computed by j www.jsoftware.com'
1-&:(1&o.d._1)0
0.459698
translated from c
)
mp=: +/ .* NB. matrix product
NB. Evaluates Simpson's Rule, also returning m and f(m) to reuse.
uquad_simpsons_mem=: adverb define
'a fa b fb'=. y
em=. a ([ + [: -: -~) b
fm=. u em
simp=. ((| b - a) % 6) * 1 4 1 mp fa , fm , fb
em, fm, simp
)
Simp=: 1 :'2{m'
Fm=: 1 :'1{m'
M=: 1 :'0{m'
NB. Efficient recursive implementation of adaptive Simpson's rule.
NB. Function values at the start, middle, end of the intervals are retained.
uquad_asr=: adverb define
'a fa b fb eps whole em fm'=. y
lt=. u uquad_simpsons_mem(a, fa, em, fm)
rt=. u uquad_simpsons_mem(em, fm, b, fb)
delta=. lt Simp + rt Simp - whole
if. (| delta) <: eps * 15 do.
lt Simp + rt Simp + delta % 15
else.
(a, fa, em, fm, (-: eps), lt Simp, lt M, lt Fm) +&(u uquad_asr) (em, fm, b, fb, (-: eps), rt Simp, rt M, rt Fm)
end.
)
NB. Integrate u from a to b using ASR with max error of eps.
quad_asr=: adverb define
'a b eps'=. y
fa=. u a
fb=. u b
t=. u uquad_simpsons_mem a, fa, b, fb
u uquad_asr a, fa, b, fb, eps, t Simp, t M, t Fm
)
echo 'Simpson''s integration of sine from 0 to 1 = ' , ": 1&o. quad_asr 0 1 1e_9 Simpson's integration of sine from 0 to 1 = 0.459698
Java
import java.util.function.Function;
public class NumericalIntegrationAdaptiveSimpsons {
public static void main(String[] args) {
Function<Double,Double> f = x -> sin(x);
System.out.printf("integrate sin(x), x = 0 .. Pi = %2.12f. Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, Math.PI, 1e-8), functionCount);
functionCount = 0;
System.out.printf("integrate sin(x), x = 0 .. 1 = %2.12f. Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, 1, 1e-8), functionCount);
}
private static double quadratureAdaptiveSimpsons(Function<Double,Double> function, double a, double b, double error) {
double fa = function.apply(a);
double fb = function.apply(b);
Triple t = quadratureAdaptiveSimpsonsOne(function, a, fa, b ,fb);
return quadratureAdaptiveSimpsonsRecursive(function, a, fa, b, fb, error, t.s, t.x, t.fx);
}
private static double quadratureAdaptiveSimpsonsRecursive(Function<Double,Double> function, double a, double fa, double b, double fb, double error, double whole, double m, double fm) {
Triple left = quadratureAdaptiveSimpsonsOne(function, a, fa, m, fm);
Triple right = quadratureAdaptiveSimpsonsOne(function, m, fm, b, fb);
double delta = left.s + right.s - whole;
if ( Math.abs(delta) <= 15*error ) {
return left.s + right.s + delta / 15;
}
return quadratureAdaptiveSimpsonsRecursive(function, a, fa, m, fm, error/2, left.s, left.x, left.fx) +
quadratureAdaptiveSimpsonsRecursive(function, m, fm, b, fb, error/2, right.s, right.x, right.fx);
}
private static Triple quadratureAdaptiveSimpsonsOne(Function<Double,Double> function, double a, double fa, double b, double fb) {
double m = (a + b) / 2;
double fm = function.apply(m);
return new Triple(m, fm, Math.abs(b-a) / 6 * (fa + 4*fm + fb));
}
private static class Triple {
double x, fx, s;
private Triple(double m, double fm, double s) {
this.x = m;
this.fx = fm;
this.s = s;
}
}
private static int functionCount = 0;
private static double sin(double x) {
functionCount++;
return Math.sin(x);
}
}
- Output:
integrate sin(x), x = 0 .. Pi = 1.999999999998. Function calls = 121 integrate sin(x), x = 0 .. 1 = 0.459697694131. Function calls = 33
Julia
Originally from Modesto Mas, https://mmas.github.io/simpson-integration-julia
function simps(f::Function, a::Number, b::Number, n::Number)
iseven(n) || throw("n must be even, and $n was given")
h = (b-a)/n
s = f(a) + f(b)
s += 4 * sum(f.(a .+ collect(1:2:n) .* h))
s += 2 * sum(f.(a .+ collect(2:2:n-1) .* h))
h/3 * s
end
println("Simpson's rule integration of sin from 0 to 1 is: ", simps(sin, 0.0, 1.0, 100))
- Output:
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936
Kotlin
// Version 1.2.71
import kotlin.math.abs
import kotlin.math.sin
typealias F = (Double) -> Double
typealias T = Triple<Double, Double, Double>
/* "structured" adaptive version, translated from Racket */
fun quadSimpsonsMem(f: F, a: Double, fa: Double, b: Double, fb: Double): T {
// Evaluates Simpson's Rule, also returning m and f(m) to reuse
val m = (a + b) / 2
val fm = f(m)
val simp = abs(b - a) / 6 * (fa + 4 * fm + fb)
return T(m, fm, simp)
}
fun quadAsrRec(f: F, a: Double, fa: Double, b: Double, fb: Double,
eps: Double, whole: Double, m: Double, fm: Double): Double {
// Efficient recursive implementation of adaptive Simpson's rule.
// Function values at the start, middle, end of the intervals are retained.
val (lm, flm, left) = quadSimpsonsMem(f, a, fa, m, fm)
val (rm, frm, right) = quadSimpsonsMem(f, m, fm, b, fb)
val delta = left + right - whole
if (abs(delta) <= eps * 15) return left + right + delta / 15
return quadAsrRec(f, a, fa, m, fm, eps / 2, left, lm, flm) +
quadAsrRec(f, m, fm, b, fb, eps / 2, right, rm, frm)
}
fun quadAsr(f: F, a: Double, b: Double, eps: Double): Double {
// Integrate f from a to b using ASR with max error of eps.
val fa = f(a)
val fb = f(b)
val (m, fm, whole) = quadSimpsonsMem(f, a, fa, b, fb)
return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}
fun main(args: Array<String>) {
val a = 0.0
val b = 1.0
val sinx = quadAsr(::sin, a, b, 1.0e-09)
println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}")
}
- Output:
Simpson's integration of sine from 0.0 to 1.0 = 0.459698
Lambdatalk
1) using the sin's primitive -cos(x)
∫01 sin(x) = {- {cos 0} {cos 1}}
-> 0.45969769413186023 // will be the reference value
2) using the adaptative simpson method
{def adasimp
{def adasimp.calc
{lambda {:fun :a :fa :b :fb}
{let { {:delta {- :b :a}}
{:m {/ {+ :a :b} 2}}
{:fa :fa} {:fb :fb}
{:fm {:fun {/ {+ :a :b} 2}}}
} {A.new :m :fm {/ {* {abs :delta} {+ :fa {* 4 :fm} :fb}} 6}} }}}
{def adasimp.rec
{lambda {:fun :a :fa :b :fb :eps :mfs}
{let { {:fun :fun}
{:a :a} {:fa :fa} {:b :b} {:fb :fb}
{:eps :eps} {:mfs :mfs}
{:left {adasimp.calc :fun :a :fa {A.get 0 :mfs} {A.get 1 :mfs}}}
{:right {adasimp.calc :fun {A.get 0 :mfs} {A.get 1 :mfs} :b :fb}}
} {let { {:fun :fun}
{:a :a} {:fa :fa} {:b :b} {:fb :fb}
{:eps :eps} {:mfs :mfs}
{:left :left} {:right :right}
{:delta {+ {A.get 2 :left}
{A.get 2 :right}
{- {A.get 2 :mfs}}}}
} {if {<= {abs :delta} {* :eps 15}}
then {+ {A.get 2 :left}
{A.get 2 :right}
{/ :delta 15}}
else {+ {adasimp.rec :fun :a :fa
{A.get 0 :mfs} {A.get 1 :mfs}
{/ :eps 2} :left }
{adasimp.rec :fun
{A.get 0 :mfs} {A.get 1 :mfs}
:b :fb {/ :eps 2} :right }}} }}}}
{lambda {:fun :a :b :eps}
{adasimp.rec :fun :a {:fun :a} :b {:fun :b} :eps
{adasimp.calc :fun :a {:fun :a} :b {:fun :b}}} }}
-> adasimp
{adasimp sin 0 1 1.0e-11}
-> 0.45969769413186023
(0.45969769413186023) // reference value
Nim
Direct translation of Python code from Wikipedia entry.
import math, sugar
type Func = float -> float
func quadSimpsonsMem(f: Func; a, fa, b, fb: float): tuple[m, fm, val: float] =
## Evaluates the Simpson's Rule, also returning m and f(m) to reuse
result.m = (a + b) / 2
result.fm = f(result.m)
result.val = abs(b - a) * (fa + 4 * result.fm + fb) / 6
func quadAsr(f: Func; a, fa, b, fb, eps, whole, m, fm: float): float =
## Efficient recursive implementation of adaptive Simpson's rule.
## Function values at the start, middle, end of the intervals are retained.
let (lm, flm, left) = f.quadSimpsonsMem(a, fa, m, fm)
let (rm, frm, right) = f.quadSimpsonsMem(m, fm, b, fb)
let delta = left + right - whole
result = if abs(delta) <= 15 * eps:
left + right + delta / 15
else:
f.quadAsr(a, fa, m, fm, eps / 2, left, lm, flm) +
f.quadAsr(m, fm, b, fb, eps / 2, right, rm, frm)
func quadAsr(f: Func; a, b, eps: float): float =
## Integrate f from a to b using Adaptive Simpson's Rule with max error of eps.
let fa = f(a)
let fb = f(b)
let (m, fm, whole) = f.quadSimpsonsMem(a, fa, b, fb)
result = f.quadAsr(a, fa, b, fb, eps, whole, m, fm)
echo "Simpson's integration of sine from 0 to 1 = ", sin.quadAsr(0, 1, 1e-9)
- Output:
Simpson's integration of sine from 0 to 1 = 0.4596976941317859
Perl
use strict;
use warnings;
sub adaptive_Simpson_quadrature {
my($f, $left, $right, $eps) = @_;
my $lf = eval "$f($left)";
my $rf = eval "$f($right)";
my ($mid, $midf, $whole) = Simpson_quadrature_mid($f, $left, $lf, $right, $rf);
return recursive_Simpsons_asr($f, $left, $lf, $right, $rf, $eps, $whole, $mid, $midf);
sub Simpson_quadrature_mid {
my($g, $l, $lf, $r, $rf) = @_;
my $mid = ($l + $r) / 2;
my $midf = eval "$g($mid)";
($mid, $midf, abs($r - $l) / 6 * ($lf + 4 * $midf + $rf))
}
sub recursive_Simpsons_asr {
my($h, $a, $fa, $b, $fb, $eps, $whole, $m, $fm) = @_;
my ($lm, $flm, $left) = Simpson_quadrature_mid($h, $a, $fa, $m, $fm);
my ($rm, $frm, $right) = Simpson_quadrature_mid($h, $m, $fm, $b, $fb);
my $delta = $left + $right - $whole;
abs($delta) <= 15 * $eps
? $left + $right + $delta / 15
: recursive_Simpsons_asr($h, $a, $fa, $m, $fm, $eps/2, $left, $lm, $flm) +
recursive_Simpsons_asr($h, $m, $fm, $b, $fb, $eps/2, $right, $rm, $frm)
}
}
my ($a, $b) = (0, 1);
my $sin = adaptive_Simpson_quadrature('sin', $a, $b, 1e-9);
printf "Simpson's integration of sine from $a to $b = %.9f", $sin
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
Phix
with javascript_semantics function quadSimpsonsMem(integer f, atom a, fa, b, fb) -- Evaluates Simpson's Rule, also returning m and f(m) to reuse. atom m = (a + b) / 2, fm = f(m), simp = abs(b-a) / 6 * (fa + 4*fm + fb) return {m, fm, simp} end function function quadAsrRec(integer f, atom a, fa, b, fb, eps, whole, m, fm) -- Efficient recursive implementation of adaptive Simpson's rule. -- Function values at the start, middle, end of the intervals are retained. atom {lm, flm, left} := quadSimpsonsMem(f, a, fa, m, fm), {rm, frm, rght} := quadSimpsonsMem(f, m, fm, b, fb), delta := left + rght - whole if abs(delta) <= eps*15 then return left + rght + delta/15 end if return quadAsrRec(f, a, fa, m, fm, eps/2, left, lm, flm) + quadAsrRec(f, m, fm, b, fb, eps/2, rght, rm, frm) end function function quadAsr(integer f, atom a, b, eps) -- Integrate f from a to b using ASR with max error of eps. atom fa := f(a), fb := f(b), {m, fm, whole} := quadSimpsonsMem(f, a, fa, b, fb) return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm) end function -- we need a mini wrapper to get a routine_id for sin() -- (because sin() is implemented in low-level assembly) function _sin(atom a) return sin(a) end function atom a := 0.0, b := 1.0, sinx := quadAsr(_sin, a, b, 1e-09) printf(1,"Simpson's integration of sine from %g to %g = %f\n", {a, b, sinx})
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Python
#! python3
'''
example
$ python3 /tmp/integrate.py
Simpson's integration of sine from 0.0 to 1.0 = 0.4596976941317858
expected answer computed by j www.jsoftware.com
1-&:(1&o.d._1)0
0.459698
translated from c
'''
import math
import collections
triple = collections.namedtuple('triple', 'm fm simp')
def _quad_simpsons_mem(f: callable, a: float , fa: float, b: float, fb: float)->tuple:
'''
Evaluates Simpson's Rule, also returning m and f(m) to reuse.
'''
m = a + (b - a) / 2
fm = f(m)
simp = abs(b - a) / 6 * (fa + 4*fm + fb)
return triple(m, fm, simp,)
def _quad_asr(f: callable, a: float, fa: float, b: float, fb: float, eps: float, whole: float, m: float, fm: float)->float:
'''
Efficient recursive implementation of adaptive Simpson's rule.
Function values at the start, middle, end of the intervals are retained.
'''
lt = _quad_simpsons_mem(f, a, fa, m, fm)
rt = _quad_simpsons_mem(f, m, fm, b, fb)
delta = lt.simp + rt.simp - whole
return (lt.simp + rt.simp + delta/15
if (abs(delta) <= eps * 15) else
_quad_asr(f, a, fa, m, fm, eps/2, lt.simp, lt.m, lt.fm) +
_quad_asr(f, m, fm, b, fb, eps/2, rt.simp, rt.m, rt.fm)
)
def quad_asr(f: callable, a: float, b: float, eps: float)->float:
'''
Integrate f from a to b using ASR with max error of eps.
'''
fa = f(a)
fb = f(b)
t = _quad_simpsons_mem(f, a, fa, b, fb)
return _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm)
def main():
(a, b,) = (0.0, 1.0,)
sinx = quad_asr(math.sin, a, b, 1e-09);
print("Simpson's integration of sine from {} to {} = {}\n".format(a, b, sinx))
main()
Raku
(formerly Perl 6)
Fairly direct translation of the Python code.
sub adaptive-Simpson-quadrature(&f, $left, $right, \ε = 1e-9) {
my $lf = f($left);
my $rf = f($right);
my ($mid, $midf, $whole) = Simpson-quadrature-mid(&f, $left, $lf, $right, $rf);
return recursive-Simpsons-asr(&f, $left, $lf, $right, $rf, ε, $whole, $mid, $midf);
sub Simpson-quadrature-mid(&g, $l, $lf, $r, $rf){
my $mid = ($l + $r) / 2;
my $midf = g($mid);
($mid, $midf, ($r - $l).abs / 6 * ($lf + 4 * $midf + $rf))
}
sub recursive-Simpsons-asr(&h, $a, $fa, $b, $fb, $eps, $whole, $m, $fm){
my ($lm, $flm, $left) = Simpson-quadrature-mid(&h, $a, $fa, $m, $fm);
my ($rm, $frm, $right) = Simpson-quadrature-mid(&h, $m, $fm, $b, $fb);
my $delta = $left + $right - $whole;
$delta.abs <= 15 * $eps
?? $left + $right + $delta / 15
!! recursive-Simpsons-asr(&h, $a, $fa, $m, $fm, $eps/2, $left, $lm, $flm) +
recursive-Simpsons-asr(&h, $m, $fm, $b, $fb, $eps/2, $right, $rm, $frm)
}
}
my ($a, $b) = 0e0, 1e0;
my $sin = adaptive-Simpson-quadrature(&sin, $a, $b, 1e-9).round(10**-9);;
say "Simpson's integration of sine from $a to $b = $sin";
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
REXX
/*REXX program performs numerical integration using adaptive Simpson's method. */
numeric digits length( pi() ) - length(.) /*use # of digits in pi for precision. */
a= 0; b= 1; f= 'SIN' /*define values for A, B, and F. */
sinx= quadAsr('SIN',a,b,"1e" || (-digits() + 1) )
say "Simpson's integration of sine from " a ' to ' b ' = ' sinx
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: pi= 3.14159265358979323846; return pi /*pi has twenty-one decimal digits. */
r2r: return arg(1) // (pi() *2) /*normalize radians ──► a unit circle, */
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadSimp: procedure; parse arg f,a,fa,b,fb; m= (a+b) / 2; interpret 'fm=' f"(m)"
simp= abs(b-a) / 6 * (fa + 4*fm + fb); return m fm simp
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadAsr: procedure; parse arg f,a,b,eps; interpret 'fa=' f"(a)"
interpret 'fb=' f"(b)"
parse value quadSimp(f,a,fa,b,fb) with m fm whole
return quadAsrR(f,a,fa,b,fb,eps,whole,m,fm)
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadAsrR: procedure; parse arg f,a,fa,b,fb,eps,whole,m,fm; frac= digits() * 3 / 4
parse value quadSimp(f,a,fa,m,fm) with lm flm left
parse value quadSimp(f,m,fm,b,fb) with rm frm right
$= left + right - whole /*calculate delta.*/
if abs($)<=eps*frac then return left + right + $/frac
return quadAsrR(f,a,fa,m,fm,eps/2, left,lm,flm) + ,
quadAsrR(f,m,fm,b,fb,eps/2,right,rm,frm)
/*──────────────────────────────────────────────────────────────────────────────────────*/
sin: procedure; parse arg x; x= r2r(x); numeric fuzz min(5, max(1, digits() -3) )
if x=pi*.5 then return 1; if x==pi * 1.5 then return -1
if abs(x)=pi | x=0 then return 0
#= x; _= x; q= x*x
do k=2 by 2 until p=#; p= #; _= - _ * q / (k * (k+1) ); #= # + _
end /*k*/
return #
- output when using the default inputs:
Simpson's integration of sine from 0 to 1 = 0.459697694131860282602
Sidef
func adaptive_Simpson_quadrature(f, left, right, ε = 1e-9) {
func quadrature_mid(l, lf, r, rf) {
var mid = (l+r)/2
var midf = f(mid)
(mid, midf, abs(r-l)/6 * (lf + 4*midf + rf))
}
func recursive_asr(a, fa, b, fb, ε, whole, m, fm) {
var (lm, flm, left) = quadrature_mid(a, fa, m, fm)
var (rm, frm, right) = quadrature_mid(m, fm, b, fb)
var Δ = (left + right - whole)
abs(Δ) <= 15*ε
? (left + right + Δ/15)
: (__FUNC__(a, fa, m, fm, ε/2, left, lm, flm) +
__FUNC__(m, fm, b, fb, ε/2, right, rm, frm))
}
var (lf = f(left), rf = f(right))
var (mid, midf, whole) = quadrature_mid(left, lf, right, rf)
recursive_asr(left, lf, right, rf, ε, whole, mid, midf)
}
var (a = 0, b = 1)
var area = adaptive_Simpson_quadrature({ .sin }, a, b, 1e-15).round(-15)
say "Simpson's integration of sine from #{a} to #{b} ≈ #{area}"
- Output:
Simpson's integration of sine from 0 to 1 ≈ 0.45969769413186
Wren
import "/fmt" for Fmt
/* "structured" adaptive version, translated from Racket */
var quadSimpsonMem = Fn.new { |f, a, fa, b, fb|
// Evaluates Simpson's Rule, also returning m and f.call(m) to reuse.
var m = (a + b) / 2
var fm = f.call(m)
var simp = (b - a).abs / 6 * (fa + 4*fm + fb)
return [m, fm, simp]
}
var quadAsrRec // recursive
quadAsrRec = Fn.new { |f, a, fa, b, fb, eps, whole, m, fm|
// Efficient recursive implementation of adaptive Simpson's rule.
// Function values at the start, middle, end of the intervals are retained.
var r1 = quadSimpsonMem.call(f, a, fa, m, fm)
var r2 = quadSimpsonMem.call(f, m, fm, b, fb)
var lm = r1[0]
var flm = r1[1]
var left = r1[2]
var rm = r2[0]
var frm = r2[1]
var right = r2[2]
var delta = left + right - whole
if (delta.abs < eps * 15) return left + right + delta/15
return quadAsrRec.call(f, a, fa, m, fm, eps/2, left, lm, flm) +
quadAsrRec.call(f, m, fm, b, fb, eps/2, right, rm, frm)
}
var quadAsr = Fn.new { |f, a, b, eps|
// Integrate f from a to b using ASR with max error of eps.
var fa = f.call(a)
var fb = f.call(b)
var r = quadSimpsonMem.call(f, a, fa, b, fb)
var m = r[0]
var fm = r[1]
var whole = r[2]
return quadAsrRec.call(f, a, fa, b, fb, eps, whole, m, fm)
}
var a = 0
var b = 1
var sinx = quadAsr.call(Fn.new { |x| x.sin }, a, b, 1e-09)
Fmt.print("Simpson's integration of sine from $d to $d = $f", a, b, sinx)
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
zkl
# "structured" adaptive version, translated from Racket
fcn _quad_simpsons_mem(f, a,fa, b,fb){
#Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
m,fm := (a + b)/2, f(m);
return(m,fm, (b - a).abs()/6*(fa + fm*4 + fb));
}
fcn _quad_asr(f, a,fa, b,fb, eps, whole, m,fm){
# Efficient recursive implementation of adaptive Simpson's rule.
# Function values at the start, middle, end of the intervals are retained.
lm,flm,left := _quad_simpsons_mem(f, a,fa, m,fm);
rm,frm,right := _quad_simpsons_mem(f, m,fm, b,fb);
delta:=left + right - whole;
if(delta.abs() <= eps*15) return(left + right + delta/15);
_quad_asr(f, a,fa, m,fm, eps/2, left , lm,flm) +
_quad_asr(f, m,fm, b,fb, eps/2, right, rm,frm)
}
fcn quad_asr(f,a,b,eps){
#Integrate f from a to b using Adaptive Simpson's Rule with max error of eps
fa,fb := f(a),f(b);
m,fm,whole := _quad_simpsons_mem(f, a,fa, b,fb);
_quad_asr(f, a,fa, b,fb, eps,whole,m,fm);
}
sinx:=quad_asr((1.0).sin.unbind(), 0.0, 1.0, 1e-09);
println("Simpson's integration of sine from 1 to 2 = ",sinx);
- Output:
Simpson's integration of sine from 1 to 2 = 0.459698