有人可以解释为什么 scipy.integrate.quad 在积分 sin(X) 时对于同样长的范围给出不同的结果吗?
我正在尝试在程序中对任意(编码时已知)函数进行数值积分 使用数值积分方法。 我正在使用 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) 来获取弧长。 我想问的是 - 是否有更好的近似/最佳实践(?)/更快的方法来做到这一点。 如果需要更多上下文,我将根据您的意愿单独发布/在此处发布上下文。
I am trying to numerically integrate an arbitrary (known when I code) function in my program
using numerical integration methods. I am using Python 2.5.2 along with SciPy's numerical integration package. In order to get a feel for it, i decided to try integrating sin(x) and observed this behavior-
>>> 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)
I find this behavior odd because -
1. In ordinary integration, integrating over the full cycle gives zero.
2. In numerical integration, this (1) isn't necessarily the case, because you may just be
approximating the total area under the curve.
In any case, either assuming 1 is True or assuming 2 is True, I find the behavior to be inconsistent. Either both integrations (-pi to pi and 0 to 2*pi) should return 0.0 (first value in the tuple is the result and the second is the error) or return 2.257...
Can someone please explain why this is happening? Is this really an inconsistency? Can someone also tell me if I am missing something really basic about numerical methods?
In any case, in my final application, I plan to use the above method to find the arc length of a function. If someone has experience in this area, please advise me on the best policy for doing this in Python.
Edit
Note
I already have the first differential values at all points in the range stored in an array.
Current error is tolerable.
End note
I have read Wikipaedia on this. As Dimitry has pointed out, I will be integrating sqrt(1+diff(f(x), x)^2) to get the Arc Length. What I wanted to ask was - is there a better approximation/ Best practice(?) / faster way to do this. If more context is needed, I'll post it separately/ post context here, as you wish.
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(6)
quad
函数是来自旧 Fortran 库的函数。 它的工作原理是根据正在积分的函数的平坦度和斜率来判断如何处理用于数值积分的步长,以便最大限度地提高效率。 这意味着,即使分析结果相同,您也可能会从一个地区得到与下一个地区略有不同的答案。毫无疑问,两个集成都应该返回零。 返回 1/(10 万亿) 的值非常接近于零! 细微的差异是由于
quad
滚动sin
并更改其步长的方式造成的。 对于您计划的任务,quad
将满足您的所有需求。编辑:
对于你正在做的事情,我认为
quad
很好。 它速度快而且相当准确。 我的最后一句话是放心使用它,除非你发现某些东西确实出了问题。 如果它没有返回无意义的答案,那么它可能工作得很好。 不用担心。The
quad
function is a function from an old Fortran library. It works by judging by the flatness and slope of the function it is integrating how to treat the step size it uses for numerical integration in order to maximize efficiency. What this means is that you may get slightly different answers from one region to the next even if they're analytically the same.Without a doubt both integrations should return zero. Returning something that is 1/(10 trillion) is pretty close to zero! The slight differences are due to the way
quad
is rolling oversin
and changing its step sizes. For your planned task,quad
will be all you need.EDIT:
For what you're doing I think
quad
is fine. It is fast and pretty accurate. My final statement is use it with confidence unless you find something that really has gone quite awry. If it doesn't return a nonsensical answer then it is probably working just fine. No worries.我认为这可能是机器精度,因为两个答案实际上都是零。
如果您想从马口中得到答案,我会将这个问题发布在 scipy 讨论区
I think it is probably machine precision since both answers are effectively zero.
If you want an answer from the horse's mouth I would post this question on the scipy discussion board
我想说 O(10^-14) 的数字实际上为零。 你的容忍度是多少?
四边形背后的算法可能不是最好的。 您可以尝试另一种集成方法,看看是否会有所改善。 五阶龙格库塔可以是一种非常好的通用技术。
这可能只是浮点数的本质:“每个计算机科学家都知道什么应该了解浮点运算”。
I would say that a number O(10^-14) is effectively zero. What's your tolerance?
It might be that the algorithm underlying quad isn't the best. You might try another method for integration and see if that improves things. A 5th order Runge-Kutta can be a very nice general purpose technique.
It could be just the nature of floating point numbers: "What Every Computer Scientist Should Know About Floating Point Arithmetic".
这个输出对我来说似乎是正确的,因为你在这里有绝对误差估计。 在普通积分和数值积分中,sin(x) 的积分值确实应该在整个周期(任何 2*pi 长度的间隔)内具有零值,并且您的结果接近该值。
要计算弧长,您应该计算 sqrt(1+diff(f(x), x)^2) 函数的积分,其中 diff(f(x), x) 是 f(x) 的导数。 另请参阅弧长
This output seems correct to me since you have absolute error estimate here. The integral value of sin(x) is indeed should have value of zero for full period (any interval of 2*pi length) in both ordinary and numeric integration and your results is close to that value.
To evaluate arc length you should calculate integral for sqrt(1+diff(f(x), x)^2) function, where diff(f(x), x) is derivative of f(x). See also Arc length
两个答案都是相同的并且正确,即在给定容差范围内为零。
Both answers are the same and correct i.e., zero within the given tolerance.
差异来自于这样的事实:sin(x)=-sin(-x) 甚至在有限精度下也是如此。 而有限精度只能近似给出sin(x)~sin(x+2*pi)。 当然,如果四元组足够聪明,能够解决这个问题,那就太好了,但它确实无法先验地知道您给出的两个间隔的积分是否相等,或者第一个结果是更好的结果。
The difference comes from the fact that sin(x)=-sin(-x) exactly even in finite precision. Whereas finite precision only gives sin(x)~sin(x+2*pi) approximately. Sure it would be nice if quad were smart enough to figure this out, but it really has no way of knowing apriori that the integral over the two intervals you give are equivalent or that the the first is a better result.