Julia vs Mathematica : Numerical integration performance

Question

Julia 绝对是新人，有一个问题要问你。

我正在通过移植我的一些 Mathematica 和 Python 代码(主要是物理学中的科学计算等)来自学一些 Julia，然后看看是什么。
到现在为止，事情还算顺利。而且很快。到目前为止。

现在，我正在模拟一个基本的锁定放大器，它本质上采用一个(可能非常复杂的)与时间相关的信号，Uin(t) ，并产生输出 Uout(t) , 锁相在某个引用频率 fref (即突出显示 Uin(t) 的分量，该分量与引用正弦波有一定的相位关系)。 描述无关紧要，重要的是它基本上是通过计算积分来实现的(为了清楚起见，我实际上在这里省略了相位):



所以，我在 Mathematica 和 Julia 中开始并测试了这个:
我定义了一个模型 Uin(t) ，传递一些参数值，然后构建Uout(t)的数组，当时 t = 0 , 范围为 fref .

Julia :我用了 QuadGK数值积分包。

T = 0.1
f = 100.
Uin(t) = sin(2pi * f * t) + sin(2pi * 2f *t)
Uout(t, fref)::Float64 = quadgk(s -> Uin(s) * sin(2pi * fref * s), t-T, t, rtol=1E-3)[1]/T
frng = 80.:1.:120.
print(@time broadcast(Uout, 0., frng))

Mathematica

T = 0.1;
f = 100.;
Uin[t_] := Sin[2 π f t] + Sin[2 π 2 f t]
Uout[t_, fref_] := NIntegrate[Sin[2 π fref s] Uin[s], {s, t - T, t}]/T
frng = Table[i, {i, 80, 120, 1}];
Timing[Table[Uout[0, fr], {fr, frng}]]

结果:

Julia 将 i7-5xxx 笔记本电脑上的电池供电时间设置为 45 到 55 秒之间的任何时间，即 很多 ，而 Mathematica 在大约 2 秒内完成了。这种差异是可怕的，老实说，令人难以置信。我知道 Mathematica 的内核中有一些非常漂亮和精致的算法，但 Julia 就是 Julia。所以，问题是:什么给了？

P.S.:设置 f和 T如 const确实将 Julia 的时间减少到 ~8-10 秒，但是 f和 T不能是 const s 在实际程序中。除此之外，还有什么明显的我遗漏了吗？

2020 年 2 月 2 日编辑:

减慢似乎是由于算法试图在值接近零时降低精度，例如见下文:对于 fref = 95，计算需要 1 整秒(!)，而对于相邻的频率值，它会立即计算(返回的结果是 (res, error) 的元组)。似乎 quadgk 函数在非常小的值处停止):

  0.000124 seconds (320 allocations: 7.047 KiB)
fref = 94.0 (-0.08637214864144352, 9.21712218998258e-6)

  1.016830 seconds (6.67 M allocations: 139.071 MiB, 14.35% gc time)
fref = 95.0 (-6.088184966010742e-16, 1.046186419361636e-16)

  0.000124 seconds (280 allocations: 6.297 KiB)
fref = 96.0 (0.1254003757465191, 0.00010132083518769636)

注意:这与我要求生产的公差无关。此外，Mathematica 默认情况下通常会达到机器精度容差，但在接近零时会有所放缓，而 numpy/scipy 只会飞过整个过程，但产生的结果不如 Mathematica 精确(在默认设置下；对此并没有太多研究)。

Answer 1

你的问题与容错的选择有关。 1e-3 的相对误差听上去还不错，但实际上是当积分接近于零时。特别是，当 fref = 80.0 时会发生这种情况。 (以及 85、90、95，而不是 100、105 等):

julia> Uout(0.0, 80.0, f, T)
1.2104987553880609e-16

引用 quadgk 的文档字符串:

(Note that it is useful to specify a positive atol in cases where norm(I) may be zero.)

让我们尝试设置一个绝对容差，例如 1e-6，然后进行比较。首先是代码(使用来自@ARamirez 的代码):

Uin(t, f) = sin(2π * f * t) + sin(4π * f * t)

function Uout(t, fref, f , T)
    quadgk(s -> Uin(s, f) * sin(2π * fref * s), t-T, t, rtol=1e-3)[1]/T
end
function Uout_new(t, fref, f , T) # with atol
    quadgk(s -> Uin(s, f) * sin(2π * fref * s), t-T, t, rtol=1e-3, atol=1e-6)[1]/T
end

然后是基准测试(为此使用 BenchmarkTools)

using BenchmarkTools
T = 0.1
f = 100.0
freqs = 80.0:1.0:120.0

@btime Uout.(0.0, $freqs, $f, $T);
6.302 s (53344283 allocations: 1.09 GiB)

@btime Uout_new.(0.0, $freqs, $f, $T);
1.270 ms (11725 allocations: 262.08 KiB)

好的，那快了 5000 倍。这可以吗？