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.

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

Translation of: zkl

<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

zkl

Translation of: Python

<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