让我们考虑以下函数
\nfunction simpson(f, a,b, n)\n h = (b-a)/n\n h2 = 2*h\n s = f(a) + 4*sum(f.((a+h):h2:(b-h))) + 2*sum(f.((a+h2):h2:(b-h2))) + f(b)\n return h*s/3\nend\nRun Code Online (Sandbox Code Playgroud)\n其中 f 是区间 [a,b] 上的实值函数,它被细分为 n(偶数 \xe2\x89\xa5 4)个等长的子区间。
\n该表达式sf(n) = simpson(x->1/x, 0.01,1, n)计算 f(x)=1/x (-ln(0.01)) 积分的近似值作为 n 的函数;误差为 sf(n)+log(0.01)。
10:10:1000 范围内的双对数图显示了一些有趣的峰值:
\n\n(事实上,该图显示了abs(sf(n) + log(0.01))因为在位置 n=320 等处: sf(n) \xe2\x89\xa4 -log(0.01) !)
如果定义 s 的行被替换为
\ns = f(a) + 4*sum([fa+k*h) for k in 1:2:(n-1)]) + 2*sum([f(a+k*h) for k in 2:2:(n-1)]) +f(b)\nRun Code Online (Sandbox Code Playgroud)\n(相应的曲线已叠加在上图中的主曲线上。)
\n为什么峰值出现在这些精确值 30、60、320、490、640、870、980 处?
\n(在 中R,相应的范围函数seq不会产生任何尖峰。)
感谢大家的贡献。\n我的问题的所有答案都试图通过重复加法过程中浮点错误的累积来解释尖峰的出现,并提出一些避免此类问题的方法。\n所有这些都是正确且有用的,但有仍然存在一些问题:
\nrange给出错误结果的精确值,而且它们毕竟相对罕见while;这对于循环来说似乎是正确的Julia,从观察中我推断出它range确实努力补偿这些错误while循环转置Julia到 \xe2\x80\x9d 任意精度计算器\xe2\x80\x9dCalc中,结果完美无缺,即图形没有尖峰range还是LinRange真的比像我一样使用列表理解更好(在什么意义上?)[f(a+k*h) for k \xe2\x80\xa6]?我已经使用了Mikaels命题并编写了一些代码来显示两个范围中的哪一个是关闭的(1);唉,我没有看到这些统计数据中出现任何模式:
\n| n | n1 | n2 |
|---|---|---|
| 30 | 1 | 1 |
| 60 | 0 | 1 |
| 320 | 1 | 0 |
| 第490章 | 1 | 0 |
| 640 | 0 | 1 |
| 第870章 | 0 | 1 |
| 980 | 0 | 1 |
| 1060 | 1 | 0 |
| 1320 | 1 | 0 |
| 1330 | 0 | 1 |
但是,即使最初的问题仍未解决,尝试准确地表达和表达人们不理解的内容仍然是一项值得的练习。感谢您的参与。
\n该a:b:c构造意味着递增,b直到超过c。\n当增量、开始或结束为浮点数时,舍入错误可能会导致四舍五入略高于或低于预期的风险。
我们可以通过添加一些仪器来检测这一点:
\nfunction simpson(f, a,b, n)\n h = (b-a)/n\n h2 = 2*h\n first_range = (a+h):h2:(b-h)\n second_range = (a+h2):h2:(b-h2)\n n1 = length(first_range)\n n2 = length(second_range)\n println(n, ' ', n1, ' ', n2)\n @assert n == n1 + n2\n s = f(a) + 4*sum(f.(first_range)) + 2*sum(f.(second_range)) + f(b)\n return h*s/3\nend\nRun Code Online (Sandbox Code Playgroud)\n这不是一个范围错误,它们的工作方式与描述的完全一样,问题是试图依赖于精确的浮点数学,这不是问题。\n对于这些特定的数字,您碰巧在 处遇到了第一个此类舍入错误n=30。
这些范围的正确工具是LinRange(a, b, n).
function simpson(f, a,b, n)\n h = (b-a)/n\n h2 = 2*h\n first_range = LinRange(a+h, b-h, n\xc3\xb72)\n second_range = LinRange(a+h2, b-h2, n\xc3\xb72)\n n1 = length(first_range)\n n2 = length(second_range)\n #println(n, ' ', n1, ' ', n2)\n @assert n == n1 + n2 # this can't fail now\n s = f(a) + 4*sum(f.(first_range)) + 2*sum(f.(second_range)) + f(b)\n return h*s/3\nend\nRun Code Online (Sandbox Code Playgroud)\n即使在 中R,它看起来也seq有一个length.out参数,我建议这是在这里使用的唯一正确的第三个参数(这相当于LinRangeJulia 中的参数)。
编辑:根据 @DNF 评论,我仔细观察,我没有意识到range实际上有一个可选length参数(很像seq),即:
first_range = range(a+h, b-h, length=n\xc3\xb72)\nRun Code Online (Sandbox Code Playgroud)\n与(参见帮助)相比,它有更多的开销,但做了一些额外的修正,LinRange其中提高了增量的准确性。需要明确的是,这并不是导致您所看到的错误的原因;而是原因。这是由于使用浮点范围的增量引起的。我认为您不太可能看到此功能有任何差异。
编辑2:由于出现了其他问题:
\n\n\n一般原因(浮点错误)无法解释范围函数给出错误结果的精确值,而且它们毕竟相对罕见
\n
不,这就是解释。您依赖于小增量的浮点数学加起来正好达到一个特定值。有时它会成功,有时则不会。选择更多“邪恶”数字(n=17 怎么样),你会看到更多。如果您想确认的话,可以逐位检查 IEEE 754 算术。我不希望找到任何模式,因为这取决于您对 a、b 和 n 的选择。
\n您也没有看到最终值稍小的情况。也请随意添加它
\n@assert last(first_range) == b-h\n@assert last(second_range) == b-h2\nRun Code Online (Sandbox Code Playgroud)\n你会看到更多的错误。
\n\n\n人们可能会认为,这个一般原因在 while 循环内变得更糟;这对于 Julia 循环来说似乎是正确的,从观察中我推断出该范围确实努力补偿这些错误
\n
你没有 while 循环,所以我假设你指的是当你这样做时[fa+k*h) for k in 1:2:(n-1)]。
LinRange不,这基本上是构建一个固定长度的穷人。它不依赖于任何浮点数来确定确切的结束条件;它是基于1:2:(n-1)都是整数,当然不会有这些问题。
\n\n正在通过适当定义的
\nrange或LinRange比使用列表理解更好(在什么意义上?)([f(a+k*h) for k \xe2\x80\xa6])?
稍微好一点是的。LinRange将使用lerpi,根据索引线性插值值(即t = i/n和a*(1-t) + t*b.range将对浮点数进行更多修正。
您会注意到这些图中的差异吗?不。
\n\n\n(在 R 中,相应的范围函数 seq 不会产生任何尖峰。)
\n
seqR 中似乎有一个推动因素seq有一个推动因素。
> seq(0,0.6, 0.3)\n[1] 0.0000000000000000000 0.2999999999999999889 0.5999999999999999778\n> seq(0,0.59999999999, 0.3)\n[1] 0.00000000000000000000 0.29999999999999998890 0.59999999998999997697\nRun Code Online (Sandbox Code Playgroud)\n基本上它会忽略您请求的增量以压缩最后一个值。哎呀!但即便如此,也只是在一定程度上;
\n> seq(0,0.5999999999, 0.3)\n[1] 0.00000000000000000000 0.29999999999999998890\nRun Code Online (Sandbox Code Playgroud)\n所以如果你在这里进行一些比 1e-10 精度更差的计算,你会在R。再次通过用长度替换结束条件而不是依赖增量来解决。