Python 对于相等的正面积和负面积，黎曼和不会产生零

Question

我编写了一个程序来使用黎曼和来逼近积分，并使用 Python 中的 matplotlib 对其进行绘图。对于 x 轴上方和下方面积相等的函数，结果面积应该为零，但我的程序输出的是一个非常小的数字。

以下代码绘制了从 -1 到 1 的奇函数 f(x) = x^3，因此面积应该为零。我的代码将其近似为 1.68065561477562 e^-15。

这是什么原因造成的？它是 delta_x、x 或 y 中的舍入误差吗？我知道我可以将值四舍五入为零，但我想知道是否还有其他问题或方法可以解决这个问题。

我试过用 Decimal.decimal class 代替 delta_x，但我得到的数字更小。

Python代码：

import matplotlib.pyplot as plt
import numpy as np

# Approximates and graphs integral using Riemann Sum


# example function: f(x) = x^3
def f_x(x):
    return x**3

# integration range from a to b with n rectangles
a, b, n = -1, 1, 1000

# calculate delta x, list of x-values, list of y-values, and approximate area under curve
delta_x = (b - a) / n

x = np.arange(a, b+delta_x, delta_x)

y = [f_x(i) for i in x]

area = sum(y) * delta_x

# graph using matplotlib
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.bar(x, y, delta_x, alpha=.5)
plt.title('a={}, b={}, n={}'.format(a, b, n))
plt.xlabel('A = {}'.format(area))
plt.show()

Answer 1

你的代码似乎对我有效。它对 f_x 返回的 y 求和，并通过将 delta_x 添加到 arange 的第二个参数来解释只有 1000 个子区间的近似误差。不幸的是，我认为舍入误差在这里起作用。

Answer 2

是的，这只是由于浮点数不准确。

为了构建分区，np.linspace() may be more appropriate since np.arange() 可能会或可能不会包括端点，具体取决于使用非整数步长时它的舍入方式。

来自 numpy.arange() 文档：

When using a non-integer step, such as 0.1, the results will often not be consistent. It is better to use linspace for these cases.

Answer 3

你需要注意，你计算的不是原始意义上的黎曼积分。您将间隔划分为 n 个箱子，然后对 n+1 个箱子求和（这里 n = 1000 但 len(x) == 1001）。所以结果可能接近你的预期，但这肯定不是到达那里的好方法。

使用 Riemann sum 您可以将间隔划分为 n 个区间，然后对这些 n 个区间的值求和。您可以选择是计算左黎曼和、右黎曼和，还是可能取中点。

import numpy as np

def f_x(x):
    return x**3

# integration range from a to b with n rectangles
a, b, n = -1, 1, 1000

delta_x = (b - a) / float(n)

x_left = np.arange(a, b, delta_x)
x_right = np.arange(a+delta_x, b+delta_x, delta_x)
x_mid = np.arange(a+delta_x/2., b+delta_x/2., delta_x)

print len(x_left),  len(x_right), len(x_mid)          ### 1000 1000 1000


area_left = f_x(x_left).sum() * delta_x
area_right = f_x(x_right).sum() * delta_x
area_mid = f_x(x_mid).sum() * delta_x

print area_left     # -0.002
print area_right    #  0.002
print area_mid      # 1.81898940355e-15

虽然中点和已经给出了很好的结果，但对于对称函数，最好选择 n 偶数，并取左右和的平均值，

print 0.5*(area_right+area_left)   # 1.76204537072e-15

这同样接近于 0。

现在值得注意的是 numpy.arange 本身会产生一些错误。更好的选择是使用 numpy.linspace

x_left = np.linspace(a, b-delta_x, n)
x_right = np.linspace(a+delta_x, b, n)
x_mid = np.linspace(a+delta_x/2., b-delta_x/2., n)

屈服

print area_left     # -0.002
print area_right    #  0.002
print area_mid      # 8.52651282912e-17
print 0.5*(area_right+area_left)   # 5.68121938382e-17

5.68121938382e-17非常接近0，之所以不完全为0，确实是floating point inaccuracies。

著名的例子是 0.1 + 0.2 - 0.3 结果是 5.5e-17 而不是 0。这是为了表明这个简单的操作引入了与黎曼积分相同的 1e-17 阶误差。

Python 对于相等的正面积和负面积，黎曼和不会产生零

Python Riemann Sum does not yield zero for equal positive and negative areas

python

matplotlib

integral

calculus