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