Numerical integration/Adaptive Simpson's method
Lyness's (1969) 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.
Write an implementation of quadrature by adaptive Simpson’s method. Use the implementation to estimate the definite integral of sin(x) from 0 to 1. Show your output.
You can use the following pseudocode, which includes Lyness's modifications 1, 2, and 3.
; Lyness'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) * (fa + 4*fm + fb)] 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) |
11l
T Triple
Float m, fm, simp
F (m, fm, simp)
.m = m
.fm = fm
.simp = simp
F _quad_simpsons_mem(f, Float a, fa, b, fb)
‘
Evaluates Simpson's Rule, also returning m and f(m) to reuse.
’
V m = a + (b - a) / 2
V fm = f(m)
V simp = abs(b - a) / 6 * (fa + 4 * fm + fb)
R Triple(m, fm, simp)
F _quad_asr(f, Float a, fa, b, fb, eps, whole, m, fm) -> Float
‘
Efficient recursive implementation of adaptive Simpson's rule.
Function values at the start, middle, end of the intervals are retained.
’
V lt = _quad_simpsons_mem(f, a, fa, m, fm)
V rt = _quad_simpsons_mem(f, m, fm, b, fb)
V delta = lt.simp + rt.simp - whole
R (I (abs(delta) <= eps * 15) {lt.simp + rt.simp + delta / 15}
E _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))
F quad_asr(f, Float a, b, eps) -> Float
‘
Integrate f from a to b using ASR with max error of eps.
’
V fa = f(a)
V fb = f(b)
V t = _quad_simpsons_mem(f, a, fa, b, fb)
R _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm)
V (a, b) = (0.0, 1.0)
V sinx = quad_asr(x -> sin(x), a, b, 1e-09)
print(‘Simpson's integration of sine from #. to #. = #.’.format(a, b, sinx))
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
ALGOL 68
...with the helper routines nested inside the main quadrature routine.
BEGIN # Numerical integration using an Adaptive Simpson's method #
PROC quad asr = ( PROC(REAL)REAL f, REAL a, b, eps )REAL:
BEGIN
MODE TRIPLE = STRUCT( REAL m, fm, simp );
# "structured" adaptive version, translated from Racket #
PROC quad simpsons mem = ( PROC(REAL)REAL f, REAL a, fa, b, fb )TRIPLE:
BEGIN # Evaluates Simpson's Rule, also returning m and f(m) to reuse. #
REAL m = ( a + b ) / 2;
REAL fm = f( m );
REAL simp = ABS ( b - a ) / 6 * ( fa + 4*fm + fb );
( m, fm, simp )
END # quad simpsons mem # ;
PROC quad asr rec = ( PROC(REAL)REAL f, REAL a, fa, b, fb, eps, whole, m, fm )REAL:
IF # 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 );
REAL delta = simp OF lt + simp OF rt - whole;
ABS delta <= eps * 15
THEN simp OF lt + simp OF rt + delta / 15
ELSE quad asr rec( f, a, fa, m, fm, eps / 2, simp OF lt, m OF lt, fm OF lt )
+ quad asr rec( f, m, fm, b, fb, eps / 2, simp OF rt, m OF rt, fm OF rt )
FI # quad asr # ;
# Integrate f from a to b using ASR with max error of eps. #
REAL fa = f( a );
REAL fb = f( b );
TRIPLE t = quad simpsons mem( f, a, fa, b, fb );
quad asr rec( f, a, fa, b, fb, eps, simp OF t, m OF t, fm OF t )
END # quad asr # ;
REAL a = 0.0, b = 1.0;
REAL sinx = quad asr( sin, a, b, 1e-09 );
print( ( "Simpson's integration of sine from ", fixed( a, -6, 3 )
, " to ", fixed( b, -6, 3 )
, " = ", fixed( sinx, -9, 6 )
, newline
)
)
END
- Output:
Simpson's integration of sine from 0.000 to 1.000 = 0.459698
ALGOL W
which is
begin % Numerical integration using an Addaptive Simpson's method %
real procedure quad_asr ( real procedure f; real value a, b, eps ) ;
begin
record triple ( real mOf, fmOf, simpOf );
% "structured" adaptive version, translated from Racket %
reference(triple) procedure quad_simpsons_mem ( real procedure f
; real value a, fa, b, fb
) ;
begin % Evaluates Simpson's Rule, also returning m and f(m) to reuse. %
real m, fm, simp;
m := ( a + b ) / 2;
fm := f( m );
simp := abs ( b - a ) / 6 * ( fa + 4*fm + fb );
triple( m, fm, simp )
end quad_simpsons_mem ;
real procedure quad_asr_rec ( real procedure f
; real value a, fa, b, fb, eps, whole, m, fm
) ;
begin
% Efficient recursive implementation of adaptive Simpson's rule. %
% Function values at the start, middle, end of the intervals are retained. %
reference(triple) lt, rt;
real delta;
lt := quad_simpsons_mem( f, a, fa, m, fm );
rt := quad_simpsons_mem( f, m, fm, b, fb );
delta := simpOf(lt) + simpOf(rt) - whole;
if abs delta <= eps * 15
then simpOf(lt) + simpOf(rt) + delta / 15
else quad_asr_rec( f, a, fa, m, fm, eps / 2, simpOf(lt), mOf(lt), fmOf(lt) )
+ quad_asr_rec( f, m, fm, b, fb, eps / 2, simpOf(rt), mOf(rt), fmOf(rt) )
end quad_asr ;
real fa, fb;
reference(triple) t;
% 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 );
quad_asr_rec( f, a, fa, b, fb, eps, simpOf(t), mOf(t), fmOf(t) )
end quad_asr ;
real a, b, sinx;
a := 0.0; b := 1.0;
sinx := quad_asr( sin, a, b, 1'-09 );
write( s_w := 0, r_format := "A", r_w := 6, r_d := 3
, "Simpson's integration of sine from ", a, " to ", b, " = "
, r_w := 12, r_d := 8
, sinx
)
end.
- Output:
Simpson's integration of sine from 0.000 to 1.000 = 0.45969769
ATS
(* Adaptive Simpson quadrature was, I believe, one of the subjects in
my undergraduate numerical analysis coursework (in the middle
1980s). It is most likely the first recursive quadrature method to
have been published, and it is very easy to write a program that
does it. *)
#include "share/atspre_staload.hats"
(* The interface is a template function, for any of the floating-point
types. *)
extern fn {tk : tkind}
quadrature_by_adaptive_simpson
(* f is an ordinary function (compiled to a C function, as
opposed to a closure with an environment), with no effect
restrictions. I make it a function rather than a closure
so we can use ordinary C functions for "f". I could have
written "-<1>" or just "->" instead of "-<fun1>" *)
(f : g0float tk -<fun1> g0float tk,
a : g0float tk,
b : g0float tk,
tol : g0float tk, (* Tolerance. *)
depth : intGte 0) (* The maximum depth of recursion. *)
(* The type of the return value is "g0float tk", or in other words
(for purposes of this example program) any of the
floating-point types "float", "double", or "ldouble" (long
double), corresponding to the floating-point types of standard
C. A plain colon means there are no effect restrictions. Had I
written ":<1>" or ":<fun1>" instead, it would have meant the
same thing. *)
: g0float tk
(* The implementation will be for any floating-point type (although it
can be overridden, either generally or for a specific type). *)
implement {tk}
quadrature_by_adaptive_simpson (f, a, b, tol, depth) =
let
(* A shorthand. Sometimes you will see typedefs written instead as
"stadef", which is a more general notation. *)
typedef real = g0float tk
(* In ATS, you can nest functions inside functions to whatever
depth you like. *)
fn
_quad_asr_simpsons (a : real,
fa : real,
b : real,
fb : real)
(* The return type is an unboxed 3-tuple of floating point
numbers. It will be passed around as a C struct. Usually
one can leave out the "@" sign, but I like to include it. A
BOXED tuple (allocated as a struct in the heap, and passed
around as a pointer) would be written '() instead of @() *)
: @(real, real, real) =
let
(* The type of 0.5 must be casted to type "real" (= "g0float
tk"). Here I will do that rather explicitly, using a
macro. Typically it is possible to let the ATS compiler
infer which type-to-type cast is necessary. See further
below where I write "g0i2f" to do the casts from integer to
floating-point type. *)
macdef one_half = g0float2float<double_kind,tk> 0.5
val m = one_half * (a + b)
val fm = f(m)
val h = b - a
in
@(m, fm, (h / g0i2f 6) * (fa + (g0i2f 4 * fm) + fb))
end
(* _quad_asr is a recursive function, so I must write "fun"
instead of "fn". *)
fun
_quad_asr {depth : nat}
(* "depth" is a hypothetical natural number, used only
during typechecking. The notations below mean that,
on each call to _quad_asr, the value that is of type
"int depth" will have to be proven to decrease
without passing zero. If this is indeed so, then the
function is guaranteed to terminate (unless it would
run out of stack space, etc.). *)
.<depth>.
(a : real,
fa : real,
b : real,
fb : real,
tol : real,
whole : real,
m : real,
fm : real,
depth : int depth) : g0float tk =
let
(* Saying "and" here means that the computations should be
performed as if at the same time, and perhaps even at the
same time. In this case that is not necessary, although
perhaps it is good as documentation. In some instances it
is REALLY what you want. *)
val @(lm, flm, left) = _quad_asr_simpsons (a, fa, m, fm)
and @(rm, frm, right) = _quad_asr_simpsons (m, fm, b, fb)
val delta = left + right - whole
and tol_ = g0f2f 0.5 * tol
in
(* I test the value of "depth" first, so the typechecker can
know that, in later branches, the value of depth is greater
than zero, and thus can be decremented. (The ATS compiler
will otherwise refuse to let you subtract 1.) I could have
tested "depth" at the same time as the later tests, by
using "+" instead of "||" to mean OR, but why that is so is
getting too technical for this example. *)
if depth = 0 then
left + right + (delta / g0i2f 15)
else if tol_ = tol || abs (delta) <= g0i2f 15 * tol then
left + right + (delta / g0i2f 15)
else
let
val left_val = _quad_asr (a, fa, m, fm, tol_, left,
lm, flm, depth - 1)
and right_val = _quad_asr (m, fm, b, fb, tol_, right,
rm, frm, depth - 1)
in
left_val + right_val
end
end
(* Here is an alternative to writing "and". Programmers in
languages such as Python will recognize the technique: *)
val @(fa, fb) = @(f(a), f(b))
val @(m, fm, whole) = _quad_asr_simpsons (a, fa, b, fb)
in
_quad_asr (a, fa, b, fb, tol, whole, m, fm, depth)
end
(* We can compile an instance of the template, AND make that instance
available to BOTH ATS and C. The
="ext#quadrature_by_adaptive_simpson_double" tells the compiler to
use "quadrature_by_adaptive_simpson_double" as the name of the C
function, rather than some mangled version of that name. *)
extern fn
quadrature_by_adaptive_simpson_double :
(* The following means "of the same type as
'quadrature_by_adaptive_simpson' instantiated for the 'double'
floating-point type". Here it serves simply as a convenient
shorthand. *)
$d2ctype (quadrature_by_adaptive_simpson<double_kind>)
= "ext#quadrature_by_adaptive_simpson_double"
implement
quadrature_by_adaptive_simpson_double (f, a, b, tol, depth) =
quadrature_by_adaptive_simpson<double_kind> (f, a, b, tol, depth)
(* Let us call quadrature_by_adaptive_simpson_double from C rather
than from ATS, simply to demonstrate the process. *)
%{
#include <math.h>
double quadrature_by_adaptive_simpson_double (void *f,
double a, double b,
double tol, int depth);
double integral_of_sine_from_0_to_1 (double tol, int depth);
double
integral_of_sine_from_0_to_1 (double tol, int depth)
{
return
quadrature_by_adaptive_simpson_double
((void *) sin, 0, 1, tol, depth);
}
%}
(* A "main" that reads no command-line arguments. It must return an
exit status. *)
implement
main () =
let
val tol = 1e-9
val depth = 1000
(* One of the ways to make the foreign function call to C. *)
val result1 = $extfcall (double, "integral_of_sine_from_0_to_1",
tol, depth)
val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
"by call to the C interface: ", result1)
(* Another way to make a foreign function call. *)
extern fn sin : double -> double = "mac#sin"
(* Calling the same compiled function, but this time from ATS
instead of C. *)
val result2 =
quadrature_by_adaptive_simpson_double
(sin, 0.0, 1.0, tol, depth)
val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
"by call to a compiled function: ", result2)
(* Using the template function directly, with single precision
numbers. *)
extern fn sinf : float -> float = "mac#sinf"
val result3 =
quadrature_by_adaptive_simpson<float_kind>
(sinf, 0.0f, 1.0f, 1e-5f, depth)
val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
"by call to the template: ", result3)
in
0 (* Exit status. *)
end
- Output:
$ patscc -std=gnu2x -g -O2 adaptive_simpson_draft_task.dats -lm && ./a.out estimate of ∫ sin x dx from 0 to 1, by call to the C interface: 0.459698 estimate of ∫ sin x dx from 0 to 1, by call to a compiled function: 0.459698 estimate of ∫ sin x dx from 0 to 1, by call to the template: 0.459698
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;
}
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
Common Lisp
(defun %%quad-asr-simpsons (f a fa b fb)
(let* ((m (/ (+ a b) 2))
(fm (funcall f m))
(h (- b a)))
;; Common Lisp supports returning multiple values at once. There
;; is no need to return them explicitly as a data structure
;; (though that also could be done).
(values m fm (* (/ h 6)
(+ (funcall f a) (* 4 fm) (funcall f b))))))
(defun %%quad-asr (f a fa b fb tol whole m fm depth)
(multiple-value-bind (lm flm left)
(%%quad-asr-simpsons f a fa m fm)
(multiple-value-bind (rm frm right)
(%%quad-asr-simpsons f m fm b fb)
(let ((delta (- (+ left right) whole))
(tol^ (/ tol 2)))
(if (or (<= depth 0)
(= tol^ tol)
(<= (abs delta) (* 15 tol)))
(+ left right (/ delta 15))
(+ (%%quad-asr f a fa m fm tol^ left
lm flm (1- depth))
(%%quad-asr f m fm b fb tol^ right
rm frm (1- depth))))))))
(defun quad-asr (f a b tol depth)
(let ((fa (funcall f a))
(fb (funcall f b)))
(multiple-value-bind (m fm whole)
(%%quad-asr-simpsons f a fa b fb)
(%%quad-asr f a fa b fb tol whole m fm depth))))
(princ "estimated definite integral of sin(x) for x from 0 to 1: ")
(princ (quad-asr #'sin 0 1 1e-9 1000))
(terpri)
- Output:
$ clisp adaptive_simpson_task.lisp estimated definite integral of sin(x) for x from 0 to 1: 0.45969772
D
import std.math;
import std.stdio;
import std.typecons;
double
quad_asr (double function (double) f, double a, double b,
double tol, int depth)
{
auto quad_asr_simpsons_ (double a, double fa, double b, double fb)
{
auto m = 0.5 * (a + b);
auto fm = f(m);
auto h = b - a;
return tuple!("m", "fm", "quadval")
(m, fm, (h / 6) * (f(a) + 4*fm + f(b)));
}
double quad_asr_ (double a, double fa, double b, double fb,
double tol, double whole, double m, double fm,
int depth)
{
//
// Please do not ask why I do not use multiple return
// statements. Finding reasons why is left as an exercise.
//
auto retval = NaN (0);
auto left = quad_asr_simpsons_ (a, fa, m, fm);
auto right = quad_asr_simpsons_ (m, fm, b, fb);
auto delta = left.quadval + right.quadval - whole;
auto tol_ = 0.5 * tol;
if (depth <= 0 || tol_ == tol || abs (delta) <= 15 * tol)
retval = left.quadval + right.quadval + (delta / 15);
else
retval = (quad_asr_ (a, fa, m, fm, tol_, left.quadval,
left.m, left.fm, depth - 1)
+ quad_asr_ (m, fm, b, fb, tol_, right.quadval,
right.m, right.fm, depth - 1));
return retval;
}
auto fa = f(a);
auto fb = f(b);
auto whole = quad_asr_simpsons_ (a, fa, b, fb);
return quad_asr_ (a, fa, b, fb, tol, whole.quadval,
whole.m, whole.fm, depth);
}
int
main ()
{
auto result = quad_asr (&sin, 0.0, 1.0, 1e-9, 1000);
printf ("estimate of ∫ sin x dx from 0 to 1: %lf\n", result);
return 0;
}
- Output:
$ gdc -g adaptive_simpson_task.d && ./a.out estimate of ∫ sin x dx from 0 to 1: 0.459698
Factor
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
- Output:
Simpson's rule integration of sin from 0 to 1 is: 0.4596976941573994
Fortran
program adaptive_simpson_task
implicit none
integer, parameter :: dbl = kind (1.0d0)
interface
function function_dbl_to_dbl (x) result (y)
import dbl
real(dbl), intent(in) :: x
real(dbl) :: y
end function function_dbl_to_dbl
end interface
write (*, 1000) quad_asr (sine, 0.0_dbl, 1.0_dbl, 1E-9_dbl, 1000)
1000 format ("estimated definite integral of sin(x) ", &
& "for x from 0 to 1: ", F15.13)
contains
function sine (x) result (y)
real(dbl), intent(in) :: x
real(dbl) :: y
y = sin (x)
end function sine
function quad_asr (f, a, b, tol, depth) result (quadval)
procedure(function_dbl_to_dbl) :: f
real(dbl), intent(in) :: a, b, tol
integer, intent(in) :: depth
real(dbl) :: quadval
real(dbl) :: fa, fb, m, fm, whole
fa = f(a)
fb = f(b)
call quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, whole)
quadval = quad_asr_ (f, a, fa, b, fb, tol, whole, m, fm, depth)
end function quad_asr
recursive function quad_asr_ (f, a, fa, b, fb, tol, whole, &
& m, fm, depth) result (quadval)
procedure(function_dbl_to_dbl) :: f
real(dbl), intent(in) :: a, fa, b, fb, tol, whole, m, fm
integer, intent(in) :: depth
real(dbl) :: quadval
real(dbl) :: lm, flm, left
real(dbl) :: rm, frm, right
real(dbl) :: delta, tol_
call quad_asr_simpsons_ (f, a, fa, m, fm, lm, flm, left)
call quad_asr_simpsons_ (f, m, fm, b, fb, rm, frm, right)
delta = left + right - whole
tol_ = tol / 2
if (depth <= 0 .or. tol_ == tol .or. abs (delta) <= 15 * tol) then
quadval = left + right + (delta / 15)
else
quadval = quad_asr_ (f, a, fa, m, fm, tol_, left, &
& lm, flm, depth - 1)
quadval = quadval + quad_asr_ (f, m, fm, b, fb, tol_, &
& right, rm, frm, depth - 1)
end if
end function quad_asr_
subroutine quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, quadval)
procedure(function_dbl_to_dbl) :: f
real(dbl), intent(in) :: a, fa, b, fb
real(dbl), intent(out) :: m, fm, quadval
m = (a + b) / 2
fm = f(m)
quadval = ((b - a) / 6) * (fa + (4 * fm) + fb)
end subroutine quad_asr_simpsons_
end program adaptive_simpson_task
- Output:
$ gfortran -O -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318
On x86-64, if you leave out the optimizer you get a trampoline:
$ gfortran -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out adaptive_simpson_task.f90:20:2: 20 | function sine (x) result (y) | ^ Warning: trampoline generated for nested function ‘sine’ [-Wtrampolines] /usr/lib/gcc/x86_64-pc-linux-gnu/12/../../../../x86_64-pc-linux-gnu/bin/ld: warning: /tmp/cc2mQOD2.o: requires executable stack (because the .note.GNU-stack section is executable) estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318
Such a bit of executable code on the stack often is not allowed on systems hardened against hackers. Otherwise you probably needn't be concerned, even though warnings are being generated both by the compiler and the linker. :)
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)
}
- Output:
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
Java
import java.util.function.Function;
public class NumericalIntegrationAdaptiveSimpsons {
public static void main(String[] args) {
Function<Double,Double> f = x -> sin(x);
System.out.printf("integrate sin(x), x = 0 .. Pi = %2.12f. Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, Math.PI, 1e-8), functionCount);
functionCount = 0;
System.out.printf("integrate sin(x), x = 0 .. 1 = %2.12f. Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, 1, 1e-8), functionCount);
}
private static double quadratureAdaptiveSimpsons(Function<Double,Double> function, double a, double b, double error) {
double fa = function.apply(a);
double fb = function.apply(b);
Triple t = quadratureAdaptiveSimpsonsOne(function, a, fa, b ,fb);
return quadratureAdaptiveSimpsonsRecursive(function, a, fa, b, fb, error, t.s, t.x, t.fx);
}
private static double quadratureAdaptiveSimpsonsRecursive(Function<Double,Double> function, double a, double fa, double b, double fb, double error, double whole, double m, double fm) {
Triple left = quadratureAdaptiveSimpsonsOne(function, a, fa, m, fm);
Triple right = quadratureAdaptiveSimpsonsOne(function, m, fm, b, fb);
double delta = left.s + right.s - whole;
if ( Math.abs(delta) <= 15*error ) {
return left.s + right.s + delta / 15;
}
return quadratureAdaptiveSimpsonsRecursive(function, a, fa, m, fm, error/2, left.s, left.x, left.fx) +
quadratureAdaptiveSimpsonsRecursive(function, m, fm, b, fb, error/2, right.s, right.x, right.fx);
}
private static Triple quadratureAdaptiveSimpsonsOne(Function<Double,Double> function, double a, double fa, double b, double fb) {
double m = (a + b) / 2;
double fm = function.apply(m);
return new Triple(m, fm, Math.abs(b-a) / 6 * (fa + 4*fm + fb));
}
private static class Triple {
double x, fx, s;
private Triple(double m, double fm, double s) {
this.x = m;
this.fx = fm;
this.s = s;
}
}
private static int functionCount = 0;
private static double sin(double x) {
functionCount++;
return Math.sin(x);
}
}
- Output:
integrate sin(x), x = 0 .. Pi = 1.999999999998. Function calls = 121 integrate sin(x), x = 0 .. 1 = 0.459697694131. Function calls = 33
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))
- Output:
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936
Kotlin
// 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
Lambdatalk
1) using the sin's primitive -cos(x)
∫01 sin(x) = {- {cos 0} {cos 1}}
-> 0.45969769413186023 // will be the reference value
2) using the adaptative simpson method
{def adasimp
{def adasimp.calc
{lambda {:fun :a :fa :b :fb}
{let { {:delta {- :b :a}}
{:m {/ {+ :a :b} 2}}
{:fa :fa} {:fb :fb}
{:fm {:fun {/ {+ :a :b} 2}}}
} {A.new :m :fm {/ {* {abs :delta} {+ :fa {* 4 :fm} :fb}} 6}} }}}
{def adasimp.rec
{lambda {:fun :a :fa :b :fb :eps :mfs}
{let { {:fun :fun}
{:a :a} {:fa :fa} {:b :b} {:fb :fb}
{:eps :eps} {:mfs :mfs}
{:left {adasimp.calc :fun :a :fa {A.get 0 :mfs} {A.get 1 :mfs}}}
{:right {adasimp.calc :fun {A.get 0 :mfs} {A.get 1 :mfs} :b :fb}}
} {let { {:fun :fun}
{:a :a} {:fa :fa} {:b :b} {:fb :fb}
{:eps :eps} {:mfs :mfs}
{:left :left} {:right :right}
{:delta {+ {A.get 2 :left}
{A.get 2 :right}
{- {A.get 2 :mfs}}}}
} {if {<= {abs :delta} {* :eps 15}}
then {+ {A.get 2 :left}
{A.get 2 :right}
{/ :delta 15}}
else {+ {adasimp.rec :fun :a :fa
{A.get 0 :mfs} {A.get 1 :mfs}
{/ :eps 2} :left }
{adasimp.rec :fun
{A.get 0 :mfs} {A.get 1 :mfs}
:b :fb {/ :eps 2} :right }}} }}}}
{lambda {:fun :a :b :eps}
{adasimp.rec :fun :a {:fun :a} :b {:fun :b} :eps
{adasimp.calc :fun :a {:fun :a} :b {:fun :b}}} }}
-> adasimp
{adasimp sin 0 1 1.0e-11}
-> 0.45969769413186023
(0.45969769413186023) // reference value
Nim
Direct translation of Python code from Wikipedia entry.
import math, sugar
type Func = float -> float
func quadSimpsonsMem(f: Func; a, fa, b, fb: float): tuple[m, fm, val: float] =
## Evaluates the Simpson's Rule, also returning m and f(m) to reuse
result.m = (a + b) / 2
result.fm = f(result.m)
result.val = abs(b - a) * (fa + 4 * result.fm + fb) / 6
func quadAsr(f: Func; a, fa, b, fb, eps, whole, m, fm: float): float =
## Efficient recursive implementation of adaptive Simpson's rule.
## Function values at the start, middle, end of the intervals are retained.
let (lm, flm, left) = f.quadSimpsonsMem(a, fa, m, fm)
let (rm, frm, right) = f.quadSimpsonsMem(m, fm, b, fb)
let delta = left + right - whole
result = if abs(delta) <= 15 * eps:
left + right + delta / 15
else:
f.quadAsr(a, fa, m, fm, eps / 2, left, lm, flm) +
f.quadAsr(m, fm, b, fb, eps / 2, right, rm, frm)
func quadAsr(f: Func; a, b, eps: float): float =
## Integrate f from a to b using Adaptive Simpson's Rule with max error of eps.
let fa = f(a)
let fb = f(b)
let (m, fm, whole) = f.quadSimpsonsMem(a, fa, b, fb)
result = f.quadAsr(a, fa, b, fb, eps, whole, m, fm)
echo "Simpson's integration of sine from 0 to 1 = ", sin.quadAsr(0, 1, 1e-9)
- Output:
Simpson's integration of sine from 0 to 1 = 0.4596976941317859
OCaml
let quad_asr f a b tol depth =
let _quad_asr_simpsons a fa b fb =
let m = 0.5 *. (a +. b) in
let fm = f m in
let h = b -. a in
(m, fm, (h /. 6.0) *. (f a +. (4.0 *. fm) +. f b))
in
let rec _quad_asr a fa b fb tol whole m fm depth =
let (lm, flm, left) = _quad_asr_simpsons a fa m fm in
let (rm, frm, right) = _quad_asr_simpsons m fm b fb in
let delta = left +. right -. whole in
let tol_ = 0.5 *. tol in
if depth <= 0 || tol_ = tol || Float.abs delta <= 15.0 *. tol then
left +. right +. (delta /. 15.0)
else
_quad_asr a fa m fm tol_ left lm flm (depth - 1)
+. _quad_asr m fm b fb tol_ right rm frm (depth - 1)
in
let fa = f a and fb = f b in
let (m, fm, whole) = _quad_asr_simpsons a fa b fb in
_quad_asr a fa b fb tol whole m fm depth
;;
print_string
"estimated definite integral of sin(x) for x from 0 to 1: ";
print_float (quad_asr sin 0.0 1.0 1.0e-9 1000);
print_newline ()
;;
- Output:
$ ocaml adaptive_simpson_task.ml estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132
Perl
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
- Output:
Simpson's integration of sine from 0 to 1 = 0.459697694
Phix
with javascript_semantics 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 = 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 := f(a), fb := 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(_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
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()
Raku
(formerly Perl 6)
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
REXX
/*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 #
- output when using the default inputs:
Simpson's integration of sine from 0 to 1 = 0.459697694131860282602
Scheme
(define (quad-asr f a b tol depth)
(letrec ;; Recursive let, so %%quad-asr can call itself.
((%%quad-asr-simpsons
(lambda (a fa b fb)
(let* ((m (/ (+ a b) 2))
(fm (f m))
(h (- b a)))
;; Scheme supports returning multiple values at
;; once. There is no need to return them explicitly as
;; a data structure (though that also could be done).
(values m fm (* (/ h 6) (+ (f a) (* 4 fm) (f b)))))))
(%%quad-asr
(lambda (a fa b fb tol whole m fm depth)
;; The R7RS standard specifies there must be
;; "let-values" and "let*-values" for receiving multiple
;; values. However, I do not want to assume its
;; presence. I will use the most widely supported
;; method, which is "call-with-values". (Typically
;; "let*-values" would be implemented as syntactic sugar
;; for the following.)
(call-with-values
(lambda () (%%quad-asr-simpsons a fa m fm))
(lambda (lm flm left)
(call-with-values
(lambda () (%%quad-asr-simpsons m fm b fb))
(lambda (rm frm right)
(let ((delta (- (+ left right) whole))
(tol^ (/ tol 2)))
(if (or (<= depth 0)
(= tol^ tol)
(<= (abs delta) (* 15 tol)))
(+ left right (/ delta 15))
(+ (%%quad-asr a fa m fm tol^ left
lm flm (- depth 1))
(%%quad-asr m fm b fb tol^ right
rm frm (- depth 1))))))))))))
(let ((fa (f a))
(fb (f b)))
(call-with-values (lambda () (%%quad-asr-simpsons a fa b fb))
(lambda (m fm whole)
(%%quad-asr a fa b fb tol whole m fm depth))))))
(display "estimated definite integral of sin(x) for x from 0 to 1: ")
(display (quad-asr sin 0 1 1e-9 1000))
(newline)
- Output:
Running it with Racket's R5RS implementation:
$ plt-r5rs adaptive_simpson_task.scm estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941317858
Running it with SCM:
$ scm -f adaptive_simpson_task.scm estimated definite integral of sin(x) for x from 0 to 1: 459.6976941317858e-3
Running it with Gambit Scheme:
$ gsi adaptive_simpson_task.scm estimated definite integral of sin(x) for x from 0 to 1: .4596976941317858
Sidef
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}"
- Output:
Simpson's integration of sine from 0 to 1 ≈ 0.45969769413186
Standard ML
fun quad_asr (tol, depth) (f, a, b) =
let
fun quad_asr_simpsons_ (a, fa, b, fb) =
let
val m = 0.5 * (a + b)
val fm = f (m)
val h = b - a
in
(m, fm, (h / 6.0) * (fa + (4.0 * fm) + fb))
end
fun quad_asr_ (a, fa, b, fb, tol, whole, m, fm, depth) =
let
val (lm, flm, left) = quad_asr_simpsons_ (a, fa, m, fm)
and (rm, frm, right) = quad_asr_simpsons_ (m, fm, b, fb)
val delta = left + right - whole
and tol_ = 0.5 * tol
in
if depth <= 0
orelse abs (tol_ - tol) <= 0.0 (* <-- SML silliness *)
orelse Real.abs (delta) <= 15.0 * tol then
left + right + (delta / 15.0)
else
quad_asr_ (a, fa, m, fm, tol_,
left, lm, flm, depth - 1)
+ quad_asr_ (m, fm, b, fb, tol_,
right, rm, frm, depth - 1)
end
val fa = f (a) and fb = f (b)
val (m, fm, whole) = quad_asr_simpsons_ (a, fa, b, fb)
in
quad_asr_ (a, fa, b, fb, tol, whole, m, fm, depth)
end
;
val quadrature = quad_asr (0.000000001, 1000);
print "estimated definite integral of sin(x) for x from 0 to 1: ";
print (Real.toString (quadrature (Math.sin, 0.0, 1.0)));
print "\n";
- Output:
$ mlton adaptive_simpson_task.sml && ./adaptive_simpson_task estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132
Wren
import "/fmt" for Fmt
/* "structured" adaptive version, translated from Racket */
var quadSimpsonMem = Fn.new { |f, a, fa, b, fb|
// Evaluates Simpson's Rule, also returning m and f.call(m) to reuse.
var m = (a + b) / 2
var fm = f.call(m)
var simp = (b - a).abs / 6 * (fa + 4*fm + fb)
return [m, fm, simp]
}
var quadAsrRec // recursive
quadAsrRec = Fn.new { |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.
var r1 = quadSimpsonMem.call(f, a, fa, m, fm)
var r2 = quadSimpsonMem.call(f, m, fm, b, fb)
var lm = r1[0]
var flm = r1[1]
var left = r1[2]
var rm = r2[0]
var frm = r2[1]
var right = r2[2]
var delta = left + right - whole
if (delta.abs < eps * 15) return left + right + delta/15
return quadAsrRec.call(f, a, fa, m, fm, eps/2, left, lm, flm) +
quadAsrRec.call(f, m, fm, b, fb, eps/2, right, rm, frm)
}
var quadAsr = Fn.new { |f, a, b, eps|
// Integrate f from a to b using ASR with max error of eps.
var fa = f.call(a)
var fb = f.call(b)
var r = quadSimpsonMem.call(f, a, fa, b, fb)
var m = r[0]
var fm = r[1]
var whole = r[2]
return quadAsrRec.call(f, a, fa, b, fb, eps, whole, m, fm)
}
var a = 0
var b = 1
var sinx = quadAsr.call(Fn.new { |x| x.sin }, a, b, 1e-09)
Fmt.print("Simpson's integration of sine from $d to $d = $f", a, b, sinx)
- Output:
Simpson's integration of sine from 0 to 1 = 0.459698
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);
}
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