Julia中range对象引发辛普森积分尖峰的原因探究

考虑如下函数：

function simpson(f, a,b, n)
  h = (b-a)/n
  h2 = 2*h
  s = f(a) + 4*sum(f.((a+h):h2:(b-h))) + 2*sum(f.((a+h2):h2:(b-h2))) + f(b)
  return h*s/3
end

其中f是区间[a,b]上的实值函数，该区间被划分为n个（n为≥4的偶数）等长子区间。

表达式sf(n) = simpson(x->1/x, 0.01,1, n)用于计算f(x)=1/x的积分近似值（精确值为-ln(0.01)），误差为sf(n)+log(0.01)。

图表 #

在范围10:10:1000上绘制双对数图时出现了一些明显的尖峰：

图表展示的是abs(sf(n) + log(0.01))，因为在n=320等位置，sf(n) ≤ -log(0.01)！

若将定义s的代码替换为以下内容，尖峰将消失：

s = f(a) + 4*sum([f(a+k*h) for k in 1:2:(n-1)]) + 2*sum([f(a+k*h) for k in 2:2:(n-1)]) +f(b)

（对应的曲线已叠加在上述图表的主曲线上。）

为何尖峰会精确出现在30、60、320、490、640、870、980这些数值处？

（在R语言中，对应的range函数seq不会产生此类尖峰。）

感谢所有贡献者。所有回复均尝试通过重复加法过程中的浮点误差积累来解释尖峰的出现，并提出了一些避免此类问题的方法。这些解释均正确且有用，但仍存在一些疑问：

• 通用原因（浮点误差）无法解释range函数产生错误结果的精确数值，且这些数值相对罕见
• 可推测该通用原因在while循环中会更严重；这在Julia循环中似乎成立，由此可推断range确实在努力补偿此类误差
• 但将我在Julia中编写的while循环转换为任意精度计算器Calc后，结果完美无缺，即图表中无尖峰
• 在合适定义的range或LinRange上进行广播，是否真的比我使用的列表推导（[f(a+k*h) for k …]）更好？具体在哪些方面？

我采用了Mikaels的提议，编写了代码来展示两个range中哪一个出现了偏差（偏差为1）；遗憾的是，我并未从这些统计数据中发现规律：

但即使初始问题未解决，尝试清晰阐述并呈现自己的疑惑仍是一项有价值的工作。感谢参与讨论。

内容的提问来源于stack exchange，提问作者Marc J Charpentier