user
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.

Pseudocode: Simpson's method, adaptive

|  |
| --- |
| ``` ; Lyness's Adaptive Simpson's Rule, Modifications 1, 2, 3 procedure simpson_rule(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 recursive_simpson(f, a, fa, b, fb, tol, whole, m, fm, depth)     lm, flm, left  := simpson_rule(f, a, fa, m, fm)     rm, frm, right := simpson_rule(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 recursive_simpson(f, a, fa, m, fm, tol', left , lm, flm, depth - 1) +                recursive_simpson(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 := simpson_rule(f, a, fa, b, fb)    return recursive_simpson(f, a, fa, b, fb, tol, whole, m, fm, depth)  ``` |
assistant
(formerly Perl 6)

**Works with**: Rakudo version 2018.10

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
```