使用 MATLAB 积分 0 阶 Bessel 函数

Question

我有一个简单的问题。我正在尝试使用 Matlab R2012a 评估 0 阶 Bessel 函数的不正确积分：

v = integral(@(x)(besselj(0, x), 0, Inf)

这给了我 v = 3.7573e+09。然而，这在理论上应该是 v = 1。当我想做的时候

v = integral(@(l)besselj(0,l), 0, 1000)

结果为 v = 1.0047。您能否简要解释一下，集成出了什么问题？以及如何正确集成 Bessel 类型的函数？

Answer 1

From the docs 对振荡函数进行不正确的积分：

q = integral(fun,0,Inf,'RelTol',1e-8,'AbsTol',1e-13)

在文档中的示例是

fun = @(x)x.^5.*exp(-x).*sin(x);

但我想在你的情况下尝试：

q = integral(@(x)(besselj(0, x),0,Inf,'RelTol',1e-8,'AbsTol',1e-13)

Answer 2

起初我怀疑对 Bessel 函数进行积分会产生有限的结果。然而，Mathematica/Wofram Alpha 显示结果是有限的，但它是not for the faint of heart。

但是，然后我被指向this site，它解释了如何正确执行此操作，并且积分的值应为 1。

我做了一些实验来验证他们陈述的正确性：

    F = @(z) arrayfun(@(y) quadgk(@(x)besselj(0,x), 0, y), z);

    z = 10:100:1e4;

    plot(z, F(z))

这给了：

如此明显，积分确实似乎收敛到 1。Wolfram Alpha 可耻！

（请注意，这是一种误导性的情节；尝试使用z = 10:1e4; 来做这件事，你会明白为什么。但是哦，原理是一样的）。

此图还准确地显示了您在 Matlab 中遇到的问题；积分的值就像一个在 1 附近的阻尼振荡，用于增加x。问题是，阻尼非常弱——如您所见，我的 z 需要一直到 10,000 才能生成此图，而振荡幅度仅降低了 ~0.5。

当您尝试通过弄乱“MaxIntervalCount”设置来进行不正确的积分时，您会得到以下信息：

    >> quadgk(@(x)besselj(0,x), 0, inf, 'maxintervalcount', 1e4) 

    Warning: Reached the limit on the maximum number of intervals in use.
    Approximate bound on error is 1.2e+009. The integral may not exist, or
    it may be difficult to approximate numerically.
    Increase MaxIntervalCount to 10396 to enable QUADGK to continue for 
    another iteration. 
    > In quadgk>vadapt at 317
      In quadgk at 216

无论您将MaxIntervalCount 设置多高；你会一直遇到这个错误。使用 quad、quadl 或类似的东西时也会发生类似的事情（这些是 R2012 integral 函数的基础）。

正如此警告和绘图所示，积分不适合通过在标准 MATLAB 中实现的任何求积方法进行精确近似（至少，据我所知）。

我相信，正如在物理论坛上所做的那样，正确的解析推导确实是获得结果的唯一方法，而无需求助于专门的求积方法。