Menu
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.
473
Riemann sums sometimes overestimate and other times underestimate. Riemann sums are approximations of the area under a curve, so they will almost always be slightly more than the actual area (an overestimation) or slightly less than the actual area (an underestimation).
If f is increasing over the whole interval [a, b] then a left-Riemann sum will also be a lower Riemann sum and a right-Riemann sum will be an upper Riemann sum; if f is decreasing, this correspondence is reversed.
A Riemann sum is an approximation of a region's area, obtained by adding up the areas of multiple simplified slices of the region. It is applied in calculus to formalize the method of exhaustion, used to determine the area of a region. This process yields the integral, which computes the value of the area exactly.
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
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
