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

From Rosetta Code
Content added Content deleted
m (added zkl header)
(→‎{{header|zkl}}: added code)
Line 34: Line 34:


=={{header|zkl}}==
=={{header|zkl}}==
{{trans|Python}}
<lang zkl></lang>
<lang zkl># "structured" adaptive version, translated from Racket
<lang zkl></lang>
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>
{{out}}
{{out}}
<pre>
<pre>
Simpson's integration of sine from 1 to 2 = 0.459698
</pre>
</pre>

Revision as of 22:39, 29 September 2018

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.

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)

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