从一个特殊函数积分的案例说开——记MATLAB R2014a符号积分bug

前言

本帖将通过一个典型的数值积分案例来介绍如何在 MATLAB 中有效地计算特殊函数（special functions）的数值积分。特殊函数的数值积分（尤其是无穷积分）比较常见的一个问题是：积分在满足收敛条件下积分结果为 NaN。如果我们使用的是 integral 或 quadgk，这一问题通常会以警告的形式给出: “Warning:Infinite or Not-a-Number value encountered,” 伴随该警告的恶果是，数值积分的结果为 NaN。

为了直观地说明问题，本帖将针对一个在数学、物理学、以及工程应用里较为常见的特殊函数 —— Modified Bessel Function of the First Kind（数学符号记作：，MATLAB 中对应的函数是 besseli），介绍如何克服积分结果为 NaN 的问题。本帖的方法虽然只是针对这一种特殊函数，但解决问题的思路适用于这一类数值计算问题，希望能对有兴趣的人有帮助！

最后，在分享本案例的过程中，我还会顺便指出 MATLAB R2014a 符号积分系统的一个计算错误，应该算是一个 bug。

特殊函数积分案例

先给出我要计算的积分：

首先，我们有必要先对的函数特性做一个基本了解（参考 http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html）。

可以看出：

• 
• 
• 是一个单调增函数
  

从 (2) 式可以得到，，类比，我们可以粗略的认为， 近似服从指数渐进性（严格的说，当 ，）。 所以， 是一个近似服从指数增长的函数。要保证这样一个无穷积分收敛，被积函数里必须要有一个衰减速度极快的因子，该因子衰减速度远快于的增长的速率，使得被积函数整体是一个快速衰减的函数，积分才有可能收敛。所以，（1）式的被积函数表达式里有一个快速衰减项，该项使得被积函数整体快速衰减（注：在面前，平方项的递增影响可以忽略不记）。

下面，我们来看看直接用 MATLAB integral 函数计算的结果（MATLAB 版本低于2012a的话可以用 quadgk 代替）。

1. integral(@(x) x.^2.*exp(-x.^2).*besseli(0,x),0,Inf)

复制代码

Warning:Infinite or Not-a-Number value encountered.

> Infunfun\private\integralCalc>iterateScalarValued at 349

  In funfun\private\integralCalc>vadapt at132

  In funfun\private\integralCalc at 83

  In integral at 88

ans =

   NaN

可以看到，直接计算的结果为NaN。

产生NaN的原因

产生NaN的原因其实很简单，因为按MATLAB的运算规则，一旦计算过程中出现了NaN，最终结果就为NaN。所以，这里一定是因为计算被积函数的函数值时得到了NaN。简单的说，当 x 很大时（如1000），快速递增单独计算的结果为Inf，而快速衰减项的结果为0，二者相乘，MATLAB 里会得到 Inf*0 = NaN。为了详细说明问题，我们先来计算一下被积函数当 x = 1000 的取值。直接用MATLAB数值计算

1. f=@(x)x^2.*exp(-x^2).*besseli(0,x)
   
2. f(1000)

复制代码

ans =

   NaN

数值计算得到结果为NaN。这正是因为单独的两项 besseli(0,x) 和 exp(-x^2) 分别计算时得到了 Inf 和 0，如下所示：

1. exp(-1000^2)

复制代码

ans =
     0

1. besseli(0,1000)

复制代码

ans =
   Inf

二者再相乘得到 NaN。既然在计算函数值的过程中就出现了NaN，难怪积分结果也是NaN。

我们再看看符号计算的结果：

1. vpa(sym('1000^2*exp(-1000^2)*besseli(0, 1000)'),4)

复制代码

ans =
8.195e-433857

按照双精度数值计算的精度，这个结果在MATLAB数值计算里可以认为是0。所以，数值计算中如果能将这个结果作为0处理，也就能避免积分结果为NaN了。

解决NaN的方法

这里介绍的方法适用于一类具备共同特征的特殊函数积分：被积函数里有一个无上界的快速递增项和一个快速衰减项，二者共同作用的结果使得被积函数在积分区间上或者积分变量 x 取值较大时快速衰减，从而积分收敛。

最简单的解决办法：既然问题出现在大 x 下，可以适当减小积分上限，将上限 Inf 换成一个有限的数，该数必须满足两个条件：1）必须足够大使得被积函数整体的取值很小（近似为0），2）又不能太大而导致数值计算里出现 Inf*0 = NaN。我们可以尝试将积分上限设置为100：

1. integral(@(x) x.^2.*exp(-x.^2).*besseli(0,x),0,100)

复制代码

ans =

   0.637956643215776

比较通用的解决办法： 有时候直接将被积函数里快速衰减因子和快速增长因子的乘积重新定义一下可能会更方便。针对之前的案例，我们可以考虑重新定exp(-x.^2).*besseli(0,x)，使得新的定义里 x 很大时直接返回结果0。为此，我们可以先大致估算一下besseli(0,x)在 x >700 时得到Inf，故我们可以定义：

1. function y = expBesseli0(x)
   
2. y = zeros(size(x));
   
3. id = x<700;
   
4. y(id) = exp(-x(id).^2).*besseli(0,x(id));

复制代码

然后，重新使用下列命令计算积分

1. integral(@(x) x.^2.*expBesseli0(x),0,Inf)

复制代码

ans =

   0.637956643215780

注意，上述代码里，积分上限依然是Inf，因为我们已经重定义了expBesseli0。

MATLAB/MuPad符号积分的错误

为了验证上述解决办法的合法性，我特意使用了 MuPad 符号积分来验证。出乎意料的是，MuPad 符号积分的结果完全错误！

1. syms x
   
2. int(x^2*exp(-x^2)*besseli(0,x),x,0,inf)

复制代码

ans =

-(3*exp(1/4))/8

符号积分的结果为负数，这显然是不可能的。被积函数在积分区间上是非负的，积分结果不可能是负的，这足以说明 MuPad 符号积分计算错误。但为了进一步验证，我们可以看看 Maple 的计算结果：



显然，Maple的结果是一个正数，与上述解决办法的结果相同。

最后，让我们来看看MuPad符号积分的更一般的 bug。积分式如下：

(3)式是(1)式的更一般化表述（将（3）式里 a、b 设置为 1即得（1）式）。让我们来看看 MATLAB、Maple、Mathematica 各自的计算结果吧。

MATLAB/MuPad

1. syms a b positive
   
2. int(x^2*exp(-a*x^2)*besseli(0,b*x),x,0,inf)

复制代码

ans =

-(exp(b^2/(4*a))*(b^2+ 2*a))/(8*a^(5/2))

MATLAB/MuPad结果里依然包含了负号，这显然是一个错误。

Maple



Mathematica


Maple 和 Mathematica 的结果相同，而且不含负号。

这应该算是 MATLAB 2014a 的一个 bug 了把。

学习了。关于积分时，找到NaN，ＬＺ确实提供了一个不错的解决问题的方法。
学习。学习。:)

跟我遇到的问题一致，果然高人！！！

yxxch 发表于 2014-9-4 22:17
跟我遇到的问题一致，果然高人！！！

你帖子里 NaN 产生的原因主要是因为你给 fsolve 的初值不合理，导致在你给的初值上，方程计算结果里 exp 函数得到了 Inf，导致方程结果为NaN，进而 fsolve 失效。只要将初值合理设置即可。解决方法我已经在你帖子里给出了：http://ilovematlab.cn/thread-303441-1-1.html

求教版主，y是一个超复杂的函数，里面也有meijerG函数，想求y关于snr的0到正无穷的积分，出现了一下错误是什么原因呢？
Attempt to reference field of non-structure array.

Error in sym/int (line 124)
   rSym = mupadmex('symobj::intdef',f.s,x.s,a.s,b.s,options);

Error in cs (line 37)
y=int(p,snr,0,100)

楼主，我也遇到了类似的问题，可是判断不出来是哪个部分 导致了Inf*0 = NaN，而且对被积函数赋值以后发现有个区域内的被积函数是一个inf+infi 一部分是NAN 一部分是正常的实数。怎样求解呢

https://www.ilovematlab.cn/thread-464365-1-1.html
(出处: MATLAB中文论坛)

版主您好！我在使用matlab里面的matlab function 模块的时候遇到问题：
Expected an integer value, found an non-integer variable  with value NaN .
这是我其中一个变量不是整数还是什么意思？我试图用vpa函数去限定我的变量，但是在matlab function 这个模块中不能用这个函数，还有其他的解决办法吗？

赞楼主~ 遇到同样的问题

楼主好赞

楼主大赞，您的那一篇向量化积分的帖子对我帮助很大，很感谢！但最近有个问题一直困扰本小白，求教楼主，对于尖锐峰的函数积分，为什么积分区间取得较大时，积分结果趋近于0？
以正态分布为例，f(E)~(E0, sigma), 在 E0-3sigma<E<E0+3sigma的范围内积分结果应该是0.9973，在30sigma范围内则该是无限趋近于1。但是再扩大区间，结果则明显不符合常理。代码如下：
E0 = 200; sigma = 0.05;
gaussian = @(E) 1./sigma./sqrt(2*pi).*exp(-power(E-E0, 2)./2./power(sigma, 2));
a(1) = integral(gaussian, E0-3.*sigma, E0+3.*sigma);
a(2) = integral(gaussian, E0-30.*sigma, E0+30.*sigma);
a(3) = integral(gaussian, E0-30.*sigma, 500);
a(4) = integral(gaussian, E0-30.*sigma, 1000);
a(5) = integral(gaussian, E0-30.*sigma, 5000);

结果如下所示：
a =

    0.9973    1.0000    1.0000    0.0000    0.0000。
求教楼主该问题出现原因是什么？在积分编程时又该如何避免？谢谢啦~

1./w.*exp(-w-(sqrt((X-li.*cos(phi)).^2+(Y+li.*sin(phi)).^2).*vc./al./2).^2./w./4).*(erf(vc.*Z./al./4./sqrt(w))-erf(vc.*(Z-da.*sec(I))./al./4./sqrt(w))).*exp(-(X-li.*cos(phi)).*vc./al./2);
请问这样的函数如何避免NaN,为二重积分，w和li为变量，同为一个衰减因子与增长因子之积