⚠️ Warning: This is a draft ⚠️
This means it might contain formatting issues, incorrect code, conceptual problems, or other severe issues.
If you want to help to improve and eventually enable this page, please fork RosettaGit's repository and open a merge request on GitHub.
{{draft task|Arithmetic operations}} Lychee (1969)'s Modified [[wp:Adaptive_Simpson's_method|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.
{| class="mw-collapsible"
| + Pseudocode: Simpson's method, adaptive |
|---|
| ''; 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) |}
C
{{trans|zkl}}
#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;
}
{{out}}
Simpson's integration of sine from 0 to 1 = 0.459698
Factor
{{trans|Julia}}
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
{{out}}
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)
}
{{out}}
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
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))
{{out}}
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936
Kotlin
{{trans|Go}}
// 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
Perl
{{trans|Perl 6}}
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
{{out}}
Simpson's integration of sine from 0 to 1 = 0.459697694
Perl 6
{{works with|Rakudo|2018.10}} 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";
{{out}}
Simpson's integration of sine from 0 to 1 = 0.459697694
Phix
{{trans|Go}}
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 = call_func(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 := call_func(f,{a}),
fb := call_func(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(routine_id("_sin"), a, b, 1e-09)
printf(1,"Simpson's integration of sine from %g to %g = %f\n", {a, b, sinx})
{{out}}
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()
REXX
{{trans|Go}}
/*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 #
{{out|output|text= when using the default inputs:}}
Simpson's integration of sine from 0 to 1 = 0.459697694131860282602
Sidef
{{trans|Perl 6}}
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}"
{{out}}
Simpson's integration of sine from 0 to 1 ≈ 0.45969769413186
zkl
{{trans|Python}}
# "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);
{{out}}
Simpson's integration of sine from 1 to 2 = 0.459698