Numerical integration/Adaptive Simpson's method: Difference between revisions

Content added Content deleted

Inline

Revision as of 01:52, 9 November 2018

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

Kotlin

Translation of: Go

<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

Works with: Rakudo version 2018.10

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

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

@@ Line 211: / Line 211: @@
 my ($a, $b) = 0e0, 1e0;
-my $sin = adaptive-Simpson-quadrature(&sin, $a, $b, 1e-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";</lang>
 {{out}}
-<pre>Simpson's integration of sine from 0 to 1 = 0.4596976941317858</pre>
+<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre>
 =={{header|zkl}}==