math - 黎曼和比高阶多项式近似收敛得更快

Question

我的函数实现的实验结果存在一些问题，我希望其他人验证我使用的函数在逻辑上是否合理。

上下文:对于某个编程/数学问题，我需要计算一个连续区间的平均值，精度为小数点后 10 位。该函数相当复杂并且涉及两个维度，因此我更愿意自己执行求和近似而不是计算连续平均值(因此必须在两个维度上对函数进行积分)。

为了帮助我，我一直在研究 Rust 中的一组逼近函数。他们计算函数 f 跨越区间 a..b 的离散平均值，使用 Riemann 和、梯形法则或固定增量 delta辛普森规则，取决于实现:

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))
}

(旁注:我在上面编写的 simpson_avg 函数实际上是 Simpson 规则的两个偏移量应用的平均值，因为它使程序不那么复杂。目前，我不相信这是问题的一部分。)

当我使 delta 接近于零时，我使用几个函数 x.atan()、x.powi(4)< 测试了每种方法的误差收敛，x，积分范围设置为 0.0 到 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

我原以为辛普森法则是所有这些函数中收敛速度最快的，但正如您所见，它的收敛速度最差，而黎曼和的表现最好。

对我来说，这毫无意义，尤其是对于简单的多项式示例，在这些示例中，辛普森规则显然可以提供更好(或至少等效)的近似值。我猜这意味着我的函数逻辑/公式存在非常微妙的问题，或者我遇到了浮点精度错误。我希望在诊断此问题时获得一些帮助。

Answer 1

正如 Shepmaster 指出的那样，您需要仔细查看迭代器走了多远。

riemann_avg 需要对 a ..= b 中的所有 x 进行迭代，然后使用两点之间的平均值(除以总和n 的 n+1 元素也是错误的)! (所以基本上 sum [ f(a+0.5*delta), f(a+1.5*delta), ..., f(b-delta/2) ])

trapezoidal_avg 只需要包括终点，否则没问题。

simpson_avg 在很多层面上都是错误的。根据wikipedia: Composite Simpson's rule你不能使用所有的三元组窗口，只能每隔一个；所以你需要奇数点，至少 3 点。

还使用了来自 https://stackoverflow.com/a/47869373/1478356 的 FloatIterator .

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