5

我正在尝试使用数值积分方法在我的程序中对任意（我编码时已知）函数进行数值积分。我正在使用 Python 2.5.2 和 SciPy 的数值集成包。为了感受它，我决定尝试整合 sin(x) 并观察这种行为 -

>>> from math import pi
>>> from scipy.integrate import quad
>>> from math import sin
>>> def integrand(x):
...     return sin(x)
... 
>>> quad(integrand, -pi, pi)
(0.0, 4.3998892617846002e-14)
>>> quad(integrand, 0, 2*pi)
(2.2579473462709165e-16, 4.3998892617846002e-14)

我觉得这种行为很奇怪，因为 -
1. 在普通积分中，整个周期的积分为零。
2. 在数值积分中，这 (1) 不一定是这种情况，因为您可能只是在近似曲线下的总面积。

无论如何，假设 1 为真或假设 2 为真，我发现行为不一致。两种积分（-pi 到 pi 和 0 到 2*pi）都应该返回 0.0（元组中的第一个值是结果，第二个是错误）或返回 2.257...

有人可以解释为什么会这样吗？这真的是矛盾吗？有人能告诉我我是否遗漏了一些关于数值方法的基本知识吗？

无论如何，在我的最终应用程序中，我打算使用上述方法来查找函数的弧长。如果有人在这方面有经验，请告诉我在 Python 中执行此操作的最佳策略。

编辑
注意
我已经在数组中存储了范围内所有点的第一个微分值。
当前误差是可以容忍的。
尾注

我已经阅读了有关此的维基百科。正如 Dimitry 所指出的，我将整合 sqrt(1+diff(f(x), x)^2) 以获得弧长。我想问的是 - 是否有更好的近似/最佳实践（？）/更快的方法来做到这一点。如果需要更多上下文，我会按照您的意愿单独发布/在此处发布上下文。

4

6 回答 6

10

该quad函数是来自旧 Fortran 库的函数。它的工作原理是通过它正在积分的函数的平坦度和斜率来判断如何处理它用于数值积分的步长，以最大限度地提高效率。这意味着您可能会从一个区域到下一个区域得到略微不同的答案，即使它们在分析上是相同的。

毫无疑问，两个积分都应该返回零。返回 1/（10 万亿）的东西非常接近于零！细微的差异是由于quad滚动方式sin和改变步长的方式。对于您计划的任务，quad将是您所需要的。

编辑：对于你正在做的事情，我认为quad这很好。它快速且非常准确。我的最后声明是自信地使用它，除非你发现某些东西确实出了问题。如果它没有返回无意义的答案，那么它可能工作得很好。不用担心。

于 2009-02-24T10:32:18.523 回答
6

我认为这可能是机器精度，因为两个答案实际上都是零。

如果您想从马口中得到答案，我会将这个问题发布在scipy 讨论板上

于 2009-02-24T10:32:55.667 回答
6

我会说一个数字 O(10^-14) 实际上为零。你的容忍度是多少？

可能是基于 quad 的算法不是最好的。您可以尝试另一种集成方法，看看是否可以改善情况。5 阶 Runge-Kutta 可能是一种非常好的通用技术。

这可能只是浮点数的本质：“每个计算机科学家都应该知道的关于浮点运算的知识”。

于 2009-02-24T10:52:04.667 回答
4

这个输出对我来说似乎是正确的，因为你在这里有绝对误差估计。在普通积分和数值积分中，sin(x) 的积分值确实应该在整个周期（2*pi 长度的任何间隔）中为零，并且您的结果接近该值。
要评估弧长，您应该计算 sqrt(1+diff(f(x), x)^2) 函数的积分，其中 diff(f(x), x) 是 f(x) 的导数。另见弧长

于 2009-02-24T10:22:12.103 回答
3 
0.0 == 2.3e-16 (absolute error tolerance 4.4e-14)

两个答案都是相同且正确的，即在给定的容差范围内为零。

于 2009-02-24T23:53:11.357 回答
2

不同之处在于 sin(x)=-sin(-x) 即使在有限精度下也是如此。而有限精度只给出了 sin(x)~sin(x+2*pi) 的近似值。当然，如果 quad 足够聪明来解决这个问题会很好，但它确实无法先验地知道你给出的两个区间上的积分是等价的，或者第一个是更好的结果。

于 2009-12-04T18:33:05.030 回答