Numerical integration/Adaptive Simpson's method

From Rosetta Code
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.

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[edit]

Translation of: 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;
}
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

Go[edit]

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

Kotlin[edit]

Translation of: 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 6[edit]

Works with: Rakudo version 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";
Output:
Simpson's integration of sine from 0 to 1 = 0.459697694

Phix[edit]

Translation of: 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})
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

zkl[edit]

Translation of: 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);
Output:
Simpson's integration of sine from 1 to 2 = 0.459698