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

Question

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

1. 首先点 P1、P2、P3、P4 的 x 和 y 位置是 计算（红点）。

2. 绿色四边形的区域是结果的一部分。 该面积是通过将其分成两个三角形（蓝色 行）。

3. 每个三角形的面积是使用 Heron 公式计算的。

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

在代码中是这样的：

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 函数只是返回二维空间中两点之间的距离。 Point 结构只是保存一个 x 和 y 坐标。 stepSize 是一个常量，我设置为 0.001

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

在互联网上我找到了使用矩形计算积分的代码：

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。

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

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

如何解释精度和步长上的这种差异？ 我的算法是不是执行错了？

更新

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

现在怎么会这样？？

Answer 1

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

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

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

这里是 Python 代码，它使用 100 个矩形以三种不同的方式逼近积分。这三种方法是trap_heron（您的方法，使用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 产生的结果几乎相同。

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

Answer 2

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