Numerical integration/Adaptive Simpson's method
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) |
C
<lang 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;
}</lang>
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Go
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article. <lang go>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)
}</lang>
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Kotlin
<lang scala>// 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)}")
}</lang>
- Output:
Simpson's integration of sine from 0.0 to 1.0 = 0.459698
Perl 6
Fairly direct translation of the Python code.
<lang perl6>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";</lang>
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
zkl
<lang 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);
}</lang> <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);</lang>
- Output:
Simpson's integration of sine from 1 to 2 = 0.459698