我试图使用数值积分方法在我的程序中数字地集成任意(当我编码时已知)函数.我使用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中执行此操作的最佳政策.
编辑
注意
我已经在存储在数组中的范围内的所有点处具有第一个差值.
当前错误是可以容忍的.
结束说明
我已经阅读了Wikipaedia.正如Dimitry指出的那样,我将整合sqrt(1 + diff(f(x),x)^ 2)来获得弧长.我想问的是 - 是否有更好的近似/最佳实践(?)/更快的方式来做到这一点.如果需要更多上下文,我会在此处单独发布/发布上下文,如您所愿.
该quad
函数是旧Fortran库中的函数.它的工作原理是通过整合功能的平坦度和斜率来整合如何处理它用于数值积分的步长,以便最大限度地提高效率.这意味着即使它们在分析上相同,您可能会从一个区域到另一个区域得到稍微不同的答案.
毫无疑问,两个集成都应该返回零.返回1 /(10万亿)的东西非常接近于零!略有不同是由于quad
滚动sin
和改变其步长的方式.对于您计划的任务,quad
将是您所需要的一切.
编辑:对于你正在做的事情,我认为没问题quad
.它快速而且非常准确.我的最后陈述是充满自信地使用它,除非你发现一些确实非常糟糕的东西.如果它没有返回荒谬的答案,那么它可能正常工作.别担心.
我认为它可能是机器精度,因为两个答案实际上都是零.
如果你想从马的口中得到答案,我会在scipy讨论板上发布这个问题
我会说数字O(10 ^ -14)实际上为零.你的宽容是什么?
可能是四边形算法不是最好的.您可以尝试其他方法进行集成,看看是否可以改进.5阶Runge-Kutta可以是一种非常好的通用技术.
它可能只是浮点数的本质:"每个计算机科学家应该知道的浮点运算".