Riemann sums are converging more quickly than higher-order polynomial approximations

Question

I'm having some issues with the experimental results of my function implementations, and I would like others to verify that the functions that I'm using are logically sound.

Context: For a certain programming/math problem, I need to calculate an average across a continuous interval to a given precision of 10 decimal places. The function is rather complicated and involves two dimensions, so I would prefer performing a sum approximation over calculating a continuous average (and therefore having to integrate the function in both dimensions) myself.

To help me out, I have been working on a set of approximation functions in Rust. They compute discrete averages of the function f across an interval a..b with a fixed increment delta using Riemann sums, the trapezoidal rule, or Simpson's rule, depending on the implementation:

pub fn riemann_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64
{
    let n = ((b - a) / delta) as usize;

    (0..n)
        .map(|i| (i as f64) * delta)
        .map(f)
        .sum::<f64>() / (n as f64)
}


pub fn trapezoidal_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64
{
    let n = ((b - a) / delta) as usize;

    (0..n)
        .map(|i| (i as f64) * delta)
        .map(f)
        .collect::<Vec<f64>>()
        .windows(2)
        .map(|xs| xs[0] + xs[1])
        .sum::<f64>() / (2.0 * (n as f64))
}


pub fn simpson_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64
{
    let n = ((b - a) / delta) as usize;

    (0..n)
        .map(|i| (i as f64) * delta)
        .map(f)
        .collect::<Vec<f64>>()
        .windows(3)
        .map(|xs| xs[0] + 4.0 * xs[1] + xs[2])
        .sum::<f64>() / (6.0 * (n as f64))
}

(Side note: The simpson_avg function I've written above is actually an average over two offset applications of Simpson's rule, since it makes the program less complex. Currently, I don't believe that this is a part of the issue.)

I tested each method's error convergence as I brought delta closer to zero, using several functions x.atan(), x.powi(4), x with the bounds of integration set from 0.0 to 1.0.

delta        riemann_err  trap_err     simpson_err
             (smaller is better)

x.atan():
0.1000000000 0.0396869230 0.0763276780 0.1136616747
0.0100000000 0.0039311575 0.0078330229 0.0117430951
0.0010000000 0.0003927407 0.0007851897 0.0011777219
0.0001000000 0.0000392703 0.0000785377 0.0001178060
0.0000100000 0.0000073928 0.0000113197 0.0000152467
0.0000010000 0.0000003927 0.0000007854 0.0000011781
0.0000001000 0.0000000393 0.0000000785 0.0000001178

x.powi(4):
0.1000000000 0.0466700000 0.0794750000 0.1081733333
0.0100000000 0.0049666670 0.0097696470 0.0145089140
0.0010000000 0.0004996667 0.0009976697 0.0014950090
0.0001000000 0.0000499967 0.0000999767 0.0001499500
0.0000100000 0.0000129997 0.0000179993 0.0000229989
0.0000010000 0.0000005000 0.0000010000 0.0000015000
0.0000001000 0.0000000500 0.0000001000 0.0000001500

x:
0.1000000000 0.0500000000 0.0950000000 0.1400000000
0.0100000000 0.0050000000 0.0099500000 0.0149000000
0.0010000000 0.0005000000 0.0009995000 0.0014990000
0.0001000000 0.0000500000 0.0000999950 0.0001499900
0.0000100000 0.0000100000 0.0000149999 0.0000199999
0.0000010000 0.0000005000 0.0000010000 0.0000015000
0.0000001000 0.0000000500 0.0000001000 0.0000001500

I expected Simpson's rule to converge the fastest out of all of these functions, but as you can see, it had the worst convergence rate, with Riemann sums performing the best.

To me, this makes no sense, especially with the simple polynomial examples where Simpson's rule clearly would provide a better (or at least equivalent) approximation. I'm guessing this means that there is either a very subtle problem with my function logic/formula, or I'm running into a floating-point precision error. I would love some help in diagnosing this problem.

If by `x` you mean the function `f(x) = x` then your output doesn't make sense. For one thing the trapezoidal error should be 0 for all delta since the area in question *is* the area of a trapezoid. Furthermore, with all of these functions, the Simpson error should be going to zero much, much faster than this. Whatever you have done with your implementation of those two methods, it doesn't work. I don't know enough Rust to say more. Why not see if you can reproduce the output of calculations that you do by hand? Floating point error doesn't explain this. You have a bug. — John Coleman, Feb 15 '18 at 01:45
FWIW, you forgot to add back `a` in all three versions. Not sure what your real problem is, though -- why don't you try to make a [mcve]? — trent, Feb 15 '18 at 02:21
@trentcl Could you elaborate on what "add back `a`" means? I've never heard of that part of that formula. Also, thank you for the suggestion; I'm working on an MCVE now. — Adam Gausmann, Feb 15 '18 at 02:48
Maybe take a look at my `FloatIterator`: https://stackoverflow.com/a/47869373/1478356 — Stefan, Feb 15 '18 at 07:09
Yeah, I didn't express that very well. Shepmaster's answer elaborates, but if you try integrating from 1.0 to 2.0 (or any interval that doesn't start at 0) it will do the wrong interval. — trent, Feb 15 '18 at 12:00
Note that you're applying Simpson's rule with an interval twice as large as either of the other approximations, so you can't expect it to converge as quickly as if you were doing the textbook version. — trent, Feb 15 '18 at 12:25

Shepmaster · Answer 1 · 2018-02-15T03:48:41.607

Think fencepost error:

fn main() {
    for i in 0..2 {
        println!("{}", i);
    }
}

0
1

You are never evaluating the last delta, therefore you always have a missing chunk.

The inclusive range operator in Rust is ..=, but it's currently unstable.

Additionally, as trentcl points out, your functions never "shift back" to account for a. This doesn't matter in your case because a is always zero, but it's still incorrect:

pub fn sample_points(a: f64, b: f64, delta: f64) {
    let n = ((b - a) / delta) as usize;

    for i in (0..n).map(|i| (i as f64) * delta) {
        println!("{}", i);
    }
}

fn main() {
    sample_points(10.0, 11.0, 0.5);
}

0
0.5

Here's the code fixing the range syntax, adding the a back into the sampling points, and avoiding the need to collect into a temporary vector:

#![feature(inclusive_range_syntax)]

extern crate itertools;
use itertools::Itertools;

pub fn riemann_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;

    (0..=n)
        .map(|i| a + (i as f64) * delta)
        .map(f)
        .sum::<f64>() / (n as f64)
}

pub fn trapezoidal_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;

    (0..=n)
        .map(|i| a + (i as f64) * delta)
        .map(f)
        .tuple_windows()
        .map(|(a, b)| a + b)
        .sum::<f64>() / (2.0 * (n as f64))
}

pub fn simpson_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;

    (0..=n)
        .map(|i| a + (i as f64) * delta)
        .map(f)
        .tuple_windows()
        .map(|(a, m, b)| a + 4.0 * m + b)
        .sum::<f64>() / (6.0 * (n as f64))
}

fn main() {
    let start = 0.;
    let end = 1.;
    let f = |x| x;
    let correct = 0.5;

    for d in 0..=7 {
        let delta = 10.0f64.powi(-d);
        let r = riemann_avg(start, end, delta, &f);
        let t = trapezoidal_avg(start, end, delta, &f);
        let s = simpson_avg(start, end, delta, &f);

        println!(
            "{:+0.8} {:+0.8} {:+0.8} {:+0.8}",
            delta,
            correct - r,
            correct - t,
            correct - s
        );
    }
}

+1.00000000 -0.50000000 +0.00000000 +0.50000000
+0.10000000 -0.05000000 -0.00000000 +0.05000000
+0.01000000 -0.00500000 +0.00000000 +0.00500000
+0.00100000 -0.00050000 -0.00000000 +0.00050000
+0.00010000 -0.00005000 -0.00000000 +0.00005000
+0.00001000 +0.00000000 +0.00000500 +0.00001000
+0.00000100 -0.00000050 -0.00000000 +0.00000050
+0.00000010 -0.00000005 +0.00000000 +0.00000005

I don't understand what's going on at `delta = 0.00001000`. My first guess is some floating point edge case, but I'm still digging. — Shepmaster, Feb 15 '18 at 03:27
In floating point arithmetic (1-0)/0.00001 = 99999.99999999999, apparently. It is better to go from a number of intervals to a delta. — red75prime, Feb 15 '18 at 06:16
`let pos = (i as f64) / (n as f64); let x = (1. - pos) * a + pos + b;` is probably more stable. — Stefan, Feb 15 '18 at 07:08
I think you've *over*corrected the Riemann sum algorithm, because you're adding up `n + 1` things but only dividing by `n`. That one really should be `0..n`. — trent, Feb 15 '18 at 12:18
@Shepmaster I've been having glitches at `delta = 0.00001` as well across several different functions, so I'm inclined to believe that may be some kind of weird floating point thing. By instead stepping by powers of two, I have had no issues like that. — Adam Gausmann, Feb 15 '18 at 17:43

Stefan · Accepted Answer · 2018-02-15T13:12:33.793

As Shepmaster pointed out you need to take a close look how far your iterator walks.

riemann_avg needs to iterator over all x in a ..= b, but then uses the average between two points (dividing the sum of n+1 elements by n would be wrong too)! (so basically sum [ f(a+0.5*delta), f(a+1.5*delta), ..., f(b-delta/2) ])

trapezoidal_avg only needed to include the end point, otherwise fine.

simpson_avg was wrong on many levels. According to wikipedia: Composite Simpson's rule you must not use all 3-tuple windows, only every second one; so you need an odd number of points, and at least 3.

Also used my FloatIterator from https://stackoverflow.com/a/47869373/1478356.

Playground

extern crate itertools;
use itertools::Itertools;

/// produces: [ linear_interpol(start, end, i/steps) | i <- 0..steps ]
///
/// linear_interpol(a, b, p) = (1 - p) * a + p * b
pub struct FloatIterator {
    current: u64,
    current_back: u64,
    steps: u64,
    start: f64,
    end: f64,
}

impl FloatIterator {
    /// results in `steps` items
    pub fn new(start: f64, end: f64, steps: u64) -> Self {
        FloatIterator {
            current: 0,
            current_back: steps,
            steps: steps,
            start: start,
            end: end,
        }
    }

    /// results in `length` items. To use the same delta as `new` increment
    /// `length` by 1.
    pub fn new_with_end(start: f64, end: f64, length: u64) -> Self {
        FloatIterator {
            current: 0,
            current_back: length,
            steps: length - 1,
            start: start,
            end: end,
        }
    }

    /// calculates number of steps from (end - start) / step
    pub fn new_with_step(start: f64, end: f64, step: f64) -> Self {
        let steps = ((end - start) / step).abs().round() as u64;
        Self::new(start, end, steps)
    }

    pub fn length(&self) -> u64 {
        self.current_back - self.current
    }

    fn at(&self, pos: u64) -> f64 {
        let f_pos = pos as f64 / self.steps as f64;
        (1. - f_pos) * self.start + f_pos * self.end
    }

    /// panics (in debug) when len doesn't fit in usize
    fn usize_len(&self) -> usize {
        let l = self.length();
        debug_assert!(l <= ::std::usize::MAX as u64);
        l as usize
    }
}

impl Iterator for FloatIterator {
    type Item = f64;

    fn next(&mut self) -> Option<Self::Item> {
        if self.current >= self.current_back {
            return None;
        }
        let result = self.at(self.current);
        self.current += 1;
        Some(result)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        let l = self.usize_len();
        (l, Some(l))
    }

    fn count(self) -> usize {
        self.usize_len()
    }
}

impl DoubleEndedIterator for FloatIterator {
    fn next_back(&mut self) -> Option<Self::Item> {
        if self.current >= self.current_back {
            return None;
        }
        self.current_back -= 1;
        let result = self.at(self.current_back);
        Some(result)
    }
}

impl ExactSizeIterator for FloatIterator {
    fn len(&self) -> usize {
        self.usize_len()
    }
}

pub fn riemann_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;
    let n = n.max(1);

    // start with:
    // [a, a+delta, ..., b-delta, b]
    // then for all neighbors (x, y) sum up f((x+y)/2)

    FloatIterator::new_with_end(a, b, n as u64 + 1)
        .tuple_windows()
        .map(|(a, b)| 0.5 * (a + b))
        .map(f)
        .sum::<f64>() / (n as f64)
}

pub fn trapezoidal_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;
    let n = n.max(1);

    // start with:
    // [a, a+delta, ..., b-delta, b]
    // then for all neighbors (x, y) sum up f((x+y)/2)

    FloatIterator::new_with_end(a, b, n as u64 + 1)
        .map(f)
        .tuple_windows()
        .map(|(a, b)| a + b)
        .sum::<f64>() / (2.0 * (n as f64))
}

pub fn simpson_avg<F>(a: f64, b: f64, delta: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let n = ((b - a) / delta) as usize;
    let n = n.max(2); // need at least 3 points in the iterator
    let n = n + (n % 2); // need odd number of points in iterator

    FloatIterator::new_with_end(a, b, n as u64 + 1)
        .map(f)
        .tuple_windows()
        .step(2)
        .map(|(a, m, b)| a + 4.0 * m + b)
        .sum::<f64>() / (3.0 * (n as f64))
}

fn compare<F, G>(a: f64, b: f64, f: F, g: G)
where
    F: Fn(f64) -> f64,
    G: Fn(f64) -> f64,
{
    let correct = g(b) - g(a);
    println!("Expected result: {:0.10}", correct);
    println!(
        "{:13} {:13} {:13} {:13}",
        "delta", "riemann_err", "trapez_err", "simpson_err"
    );
    for d in 0..8 {
        let delta = 10.0f64.powi(-d);
        let r = riemann_avg(a, b, delta, &f);
        let t = trapezoidal_avg(a, b, delta, &f);
        let s = simpson_avg(a, b, delta, &f);

        println!(
            "{:+0.10} {:+0.10} {:+0.10} {:+0.10}",
            delta,
            correct - r,
            correct - t,
            correct - s,
        );
    }
}

fn main() {
    let start = 0.;
    let end = 1.;

    println!("f(x) = atan(x)");
    compare(
        start,
        end,
        |x| x.atan(),
        |x| x * x.atan() - 0.5 * (1. + x * x).ln(),
    );

    println!("");

    println!("f(x) = x^4");
    compare(start, end, |x| x.powi(4), |x| 0.2 * x.powi(5));

    println!("");

    println!("f(x) = x");
    compare(start, end, |x| x, |x| 0.5 * x * x);
}

f(x) = atan(x)
Expected result: 0.4388245731
delta         riemann_err   trapez_err    simpson_err  
+1.0000000000 -0.0248230359 +0.0461254914 -0.0011735268
+0.1000000000 -0.0002086380 +0.0004170148 -0.0000014072
+0.0100000000 -0.0000020834 +0.0000041667 -0.0000000001
+0.0010000000 -0.0000000208 +0.0000000417 -0.0000000000
+0.0001000000 -0.0000000002 +0.0000000004 +0.0000000000
+0.0000100000 -0.0000000000 +0.0000000000 -0.0000000000
+0.0000010000 -0.0000000000 +0.0000000000 -0.0000000000
+0.0000001000 +0.0000000000 +0.0000000000 -0.0000000000

f(x) = x^4
Expected result: 0.2000000000
delta         riemann_err   trapez_err    simpson_err  
+1.0000000000 +0.1375000000 -0.3000000000 -0.0083333333
+0.1000000000 +0.0016637500 -0.0033300000 -0.0000133333
+0.0100000000 +0.0000166664 -0.0000333330 -0.0000000013
+0.0010000000 +0.0000001667 -0.0000003333 -0.0000000000
+0.0001000000 +0.0000000017 -0.0000000033 +0.0000000000
+0.0000100000 +0.0000000000 -0.0000000000 -0.0000000000
+0.0000010000 +0.0000000000 -0.0000000000 -0.0000000000
+0.0000001000 +0.0000000000 +0.0000000000 +0.0000000000

f(x) = x
Expected result: 0.5000000000
delta         riemann_err   trapez_err    simpson_err  
+1.0000000000 +0.0000000000 +0.0000000000 +0.0000000000
+0.1000000000 -0.0000000000 -0.0000000000 -0.0000000000
+0.0100000000 +0.0000000000 +0.0000000000 +0.0000000000
+0.0010000000 +0.0000000000 +0.0000000000 +0.0000000000
+0.0001000000 -0.0000000000 -0.0000000000 +0.0000000000
+0.0000100000 -0.0000000000 -0.0000000000 +0.0000000000
+0.0000010000 -0.0000000000 -0.0000000000 +0.0000000000
+0.0000001000 +0.0000000000 +0.0000000000 -0.0000000000

You don't need the fake point if you swap `.tuples()` and `.tuple_windows()`. Then your closure argument becomes `((a, m), (m, b))`. On the other hand, perhaps `.tuple_windows().step(2).map(|(a, m, b)| ...)` would be even better. — trent, Feb 15 '18 at 12:55
@trentcl right, leads to even better results (I suspect `b + delta` was the problem). — Stefan, Feb 15 '18 at 13:09

Riemann sums are converging more quickly than higher-order polynomial approximations

2 Answers2