有人可以解释为什么 scipy.integrate.quad 在积分 sin(X) 时对于同样长的范围给出不同的结果吗?

问题描述 投票:0回答:6

我正在尝试在程序中对任意(编码时已知)函数进行数值积分 使用数值积分方法。我正在使用 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 为 True 还是假设 2 为 True,我发现行为都是不一致的。两种积分(-pi 到 pi 和 0 到 2*pi)都应返回 0.0(元组中的第一个值是结果,第二个值是错误)或返回 2.257...

有人可以解释为什么会发生这种情况吗?这真的是矛盾吗?有人也可以告诉我我是否遗漏了一些关于数值方法的基本知识吗?

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

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

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

python scipy numerical-methods numerical-integration
6个回答
10
投票

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

毫无疑问,两个集成都应该返回零。返回 1/(10 万亿) 的值非常接近于零!细微的差异是由于

quad
滚动
sin
并更改其步长的方式造成的。对于您计划的任务,
quad
将满足您的所有需求。

编辑: 对于你正在做的事情,我认为

quad
很好。它速度快而且相当准确。我的最后一句话是放心使用它,除非你发现某些东西确实出了问题。如果它没有返回无意义的答案,那么它可能工作得很好。不用担心。


6
投票

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

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


6
投票

我想说,数字 O(10^-14) 实际上为零。你的容忍度是多少?

四元组底层的算法可能不是最好的。您可以尝试另一种集成方法,看看是否会有所改善。五阶龙格库塔可以是一种非常好的通用技术。

这可能只是浮点数的本质:“每个计算机科学家应该了解浮点算术”。


4
投票

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


3
投票
0.0 == 2.3e-16 (absolute error tolerance 4.4e-14)

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


2
投票

差异来自于这样的事实:sin(x)=-sin(-x) 在有限精度下完全相等。而有限精度只能近似给出sin(x)~sin(x+2*pi)。当然,如果四边形足够聪明,能够解决这个问题,那就太好了,但它确实无法先验地知道您给出的两个间隔的积分是否相等,或者第一个结果是更好的结果。

© www.soinside.com 2019 - 2024. All rights reserved.