积分计算的替代算法不准确

Question

我想到了一种积分计算算法，它应该比常规矩形方法更准确。我的算法可以用图形来最好地描述（我使用 f(x) = sin(x) 作为示例）：

首先点P1,P2,P3,P4的x和y位置为计算（红点）。
绿色四边形的面积是结果的一部分。这个面积是通过将它分成两个三角形（蓝色行).
每个三角形的面积使用海伦公式计算。

我认为这显然会比矩形方法产生更好的结果。

在代码中看起来像这样：

double integral(double f(double x), double min, double max) {
    Point p1, p2, p3, p4;
    double area = 0.0;
    double s1 = 0.0;
    double s2 = 0.0;


    for(double x = min; x < max; x += stepSize) {
        p1.x = x;
        p1.y = 0.0;

        p2.x = x;
        p2.y = f(x);

        p3.x = x + stepSize;
        p3.y = f(x + stepSize);

        p4.x = x + stepSize;
        p4.y = 0.0;


        s1 = 0.5 * (distance(p1, p2) + distance(p2, p4) + distance(p1, p4));
        s2 = 0.5 * (distance(p2, p3) + distance(p3, p4) + distance(p2, p4));

        area += sqrt(s1 * (s1 - distance(p1, p2)) * (s1 - distance(p2, p4)) * (s1 - distance(p1, p4)));
        area += sqrt(s2 * (s2 - distance(p2, p3)) * (s2 - distance(p3, p4)) * (s2 - distance(p2, p4)));
    }

    return area;
}

distance 函数只是返回二维空间中两点之间的距离 space。 Point 结构仅包含 x 和 y 坐标。 stepSize 是一个常量，我设置为 0.001

我的函数给出的结果是正确的，但我想知道它与矩形方法相比要精确多少。

在 Internet 上我找到了这段使用矩形计算积分的代码：

double integral2(double(*f)(double x), double a, double b, int n) {
    double step = (b - a) / n;  // width of each small rectangle
    double area = 0.0;  // signed area
    for (int i = 0; i < n; i ++) {
        area += f(a + (i + 0.5) * step) * step; // sum up each small rectangle
    }
    return area;
}

我都使用标准 math.h sin 函数从 0.0 到半 π 对它们进行了测试。这个面积应该是1。

我的算法给出了步长 0.001 的结果 1.000204。

矩形算法给了我结果 1.000010，计算步长为 0.015708。

如何解释精度和步长的这种差异？ 我的算法实现有误吗？

更新

使用第二种方法的计算步长，我得到的结果 0.999983 比步长 0.001 的结果更接近 1。

那怎么行呢？？

Answer 1

您可以尝试使用 Kahan 求和来减少误差，但精度问题确实存在。毕竟你是在用数值方法逼近积分。

Answer 2

您的最后一个梯形可能太宽：如果 max-min 不是 stepSize 的倍数，x+stepSize 可能会大于 max。这就是为什么在您包含的矩形求和代码中，他们使用 n（矩形的数量）而不是 stepSize。

你计算梯形的方法很复杂。注意它的面积是stepSize * (P2.y + P3.y)/2。这增加了计算成本，但我想这不是测试积分中数值错误的原因。

除了这些问题，你的方法在其他方面等同于梯形法则。 https://en.wikipedia.org/wiki/Trapezoidal_rule

这里是 Python 代码，它使用 100 个矩形以三种不同的方式逼近积分。三种方式分别是trap_heron（你的方法，使用赫伦法则）、trap（梯形法）和rect（矩形求和）。你的问题是C++，但是结果应该是一样的。

import math

N = 100

def dist(a, b):
    dx = a[0] - b[0]
    dy = a[1] - b[1]
    return math.sqrt(dx*dx + dy*dy)

def trap_heron(f, min, max):
    area = 0.0
    for i in range(N):
        x0 = min + (max-min) * i/N
        x1 = min + (max-min) * (i+1)/N
        y0 = f(x0)
        y1 = f(x1)

        p1 = (x0, 0.0)
        p2 = (x0, y0)
        p3 = (x1, y1)
        p4 = (x1, 0.0)

        s1 = 0.5 * (dist(p1, p2) + dist(p2, p4) + dist(p1, p4))
        s2 = 0.5 * (dist(p2, p3) + dist(p3, p4) + dist(p2, p4))

        area += math.sqrt(s1 * (s1 - dist(p1, p2)) * (s1 - dist(p2, p4)) * (s1 - dist(p1, p4)))
        area += math.sqrt(s2 * (s2 - dist(p2, p3)) * (s2 - dist(p3, p4)) * (s2 - dist(p2, p4)))
    return area

def trap(f, min, max):
    area = 0.0
    for i in range(N):
        x0 = min + (max-min) * i/N
        x1 = min + (max-min) * (i+1)/N
        y0 = f(x0)
        y1 = f(x1)
        area += (x1-x0) * (y0+y1)/2

    return area


def rect(f, min, max):
    area = 0.0
    for i in range(N):
        y = f(min + (max-min)*(i+0.5)/N)
        area += (max-min)/N * y
    return area

print(trap(math.sin, 0, math.pi/2))
print(trap_heron(math.sin, 0, math.pi/2))
print(rect(math.sin, 0, math.pi/2))

输出为：

0.9999794382396076
0.9999794382396054
1.0000102809119051

请注意 trap 和 trap_heron 产生的结果几乎相同。

在您的评论中，您的结果是 1.015686。误差非常接近 stepSize * sin(pi/2)，所以我猜你总结了太多的梯形。

积分计算的替代算法不准确

Alternative algorithm for integral calculation inaccurate

c++

algorithm

integral