Numerical integration/Adaptive Simpson's method
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.
; 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)
|
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