从一个特殊函数积分的案例说开——记MATLAB/MuPad R2014a又一bug
本帖最后由 winner245 于 2014-7-15 02:14 编辑前言
本帖将通过一个典型的数值积分案例来介绍如何在MATLAB中有效地计算特殊函数的数值积分。特殊函数的数值积分中（尤其是无穷积分）比较常见的一个问题是：在积分收敛的前提下积分结果为 NaN。如果我们使用的是 integral 或 quadgk 这类数值积分函数，这一问题通常会以警告的形式给出: “Warning:Infinite or Not-a-Number value encountered,” 伴随该警告的恶果是，数值积分的结果为 NaN。
为了直观地说明问题，本帖将针对一个在数学、物理学、以及工程应用里较为常见的特殊函数 -----Modified Bessel Function of the First Kind（数学符号记作：I_{\nu}(x) ，MATLAB 中对应的函数是 besseli），介绍如何有效的克服积分结果为 NaN 的问题。本帖讲述的方法虽然只是针对这一种特殊函数，但解决问题的思路适用于这一类数值计算问题，希望能对有兴趣的人有帮助！
最后，在分享这个案例的过程中，我还会顺带指出 MATLAB/MuPad 符号积分系统的一个计算错误，该错误已在 MATLAB R2014a 下验证过，应该算是一个 bug。
数值积分案例
给出我要计算的积分：\int_0^\infty x^2\mathrm{e}^{-x^2}I_0(x) \mathrm{d}x ~~~~~~~~~~~~ (1)
首先，我们有必要先对I_0(x)的函数特性做一个基本了解（参考：http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html）。
I_0(x)=\sum_{k=0}^\infty\frac{\left(\frac{x}{2}\right)^{2k}}{(k!)^2} ~~~~~~~~~~~~ (2)
可以看出：
[*]I_0(x)\geq 1 \forall x\geq 0
[*]I_0(0)=1, I_0(+\infty)=+\infty
[*]I_0(x), x\geq 0 是一个单调增函数
从 (2) 式可以得到，I_0(2x)=\sum_{k=0}^\infty \left(\frac{x^k}{k!}\right)^2，类比 \mathrm{e}^x=\sum_{k=0}^\infty \frac{x^k}{k!}，我们可以粗略的认为，I_0(2x) 近似服从指数渐进性（严格的说，当 x\to\infty，I_0(x)\sim \mathrm{e}^x/\sqrt{2\pi x}）。 所以，I_0(x) 是一个近似服从指数增长的函数。要保证这样一个无穷积分收敛，被积函数里必须要有一个衰减速度极快的因子，该因子衰减速度远快于I_0(x)的增长的速率，使得被积函数整体是一个快速衰减的函数，积分才有可能收敛。所以，（1）式的被积函数表达式里有一个快速衰减项\mathrm{e}^{-x^2}，该项使得被积函数整体快速衰减（注：在\mathrm{e}^{-x^2}面前，平方项x^2的递增影响可以忽略不记）。
下面，我们来看看直接用 integral 函数计算的结果（版本低于2012a的可以用 quadgk 代替）。 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 349In funfun\private\integralCalc>vadapt at132In funfun\private\integralCalc at 83In integral at 88 ans = NaN
可以看到，直接计算的结果为NaN。
产生NaN的原因
产生NaN的原因其实很简单，因为按MATLAB的运算规则，一旦计算过程中出现了NaN，最终结果就为NaN。所以，这里一定是因为计算被积函数的函数值时得到了NaN。简单的说，当x很大时（如1000），快速递增单独计算的结果为Inf，而快速衰减项的结果为0，二者相乘，MATLAB 里会得到Inf*0 = NaN。为了详细说明问题，我们先来计算一下被积函数当x=1000的取值。直接用MATLAB数值计算
>> f=@(x)x^2.*exp(-x^2).*besseli(0,x)>> f(1000)
ans = NaN
数值计算得到结果为NaN。这正是因为单独的两项besseli(0,x) 和exp(-x^2)分别计算时得到了Inf和0，如下所示：exp(-1000^2) ans = 0 besseli(0,1000) ans = Inf
二者再相乘得到 NaN。既然在计算函数值的过程中就出现了NaN，难怪积分结果也是NaN。
我们再看看符号计算的结果： vpa(sym('1000^2*exp(-1000^2)*besseli(0, 1000)'),4) ans = 8.195e-433857
符号计算结果为8.195e-433857，按照MATLAB双精度数值计算的精度，这个结果在数值计算里就是0。所以，如果数值计算中，能将这个结果作为0处理，也就能避免积分结果为NaN了。
解决NaN问题的方法
这里介绍的方法适用于一类具备共同特征的特殊函数积分：被积函数里有一个无上界的快速递增项和一个快速衰减项，二者共同作用的结果使得被积函数在积分区间上或者积分变量x 取值较大时快速衰减，从而积分收敛。
最简单的解决办法：既然问题出现在大x下，可以适当减小积分上限，将上限Inf换成一个有限的数，该数必须满足两个条件：1）x 必须足够大使得被积函数整体的取值很小（近似为0），2）x又不能太大而导致数值计算里出现Inf*0 = NaN。我们可以尝试将上限设置为100： 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，故我们可以定义：
function y = expBesseli0(x)y = zeros(size(x));id = x<700;y(id)= exp(-x(id).^2).*besseli(0,x(id));
然后，重新使用下列命令计算积分 integral(@(x) x.^2.*expBesseli0(x),0,Inf) ans = 0.637956643215780
注意，上述代码里，积分上限依然是Inf，因为我们已经重定义了expBesseli0。
MATLAB/MuPad符号积分的错误
为了验证上面数值计算的合法性，我特意使用了MuPad符号积分来验证。出乎意料的是，MuPad 符号积分的结果完全错误！
>> syms x>> int(x^2*exp(-x^2)*besseli(0,x),x,0,inf)
ans = -(3*exp(1/4))/8
上面的符号积分的结果居然为负数，这显然是错误的。被积函数在积分区间上是非负的，积分结果不可能是负的，这足以说明MuPad计算错误。但为了进一步验证，我又专门通过Maple来验证了一下。
最后，让我们来看看MuPad符号积分的更一般的 bug。积分式如下：
\int_0^\infty x^2\mathrm{e}^{-ax^2}I_0(bx)\mathrm{d}x, ~ a>0,~~ b>0 ~~~~~~~ (3)
(3)式是(1)式的更一般化表述（将（3）式里 a、b 设置为 1即得（1）式）。让我们来看看 MATLAB、Maple、Mathematica 各自的计算结果吧。
MATLAB/MuPad>> syms a b positive>> 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))
MuPad 结果里依然包含了负号，这显然是不可能的，因为被积函数在积分区间上非负。再让我们看看Maple和Mathematica的结果。
Maple：
Mathematica：
Maple 和 Mathematica 的结果相同，而且不含负号。
这应该算是 MATLAB 2014a 的一个 bug 了把。
建议加个精华，讨论mupad的帖子比较少，bug就更没见过了，应该够精华条件。winner老弟相当谨慎，给出如此翔实的例子论证，这几天有点事，后面验证再跟进讨论，各位继续。 bainhome 发表于 2014-7-3 00:12
建议加个精华，讨论mupad的帖子比较少，bug就更没见过了，应该够精华条件。winner老弟相当谨慎，给出如此翔 ...
呵呵，之前还发现一个 MuPad 积分的 bug：http://forum.simwe.com/thread-1106018-1-1.html。貌似是在计算 -1.3 次方时出的问题
syms x
int((-1/x)^1.3,x)
int((-x)^(-1.3),x)
ans =
-(10*(-1/x)^(3/10))/7
ans =
10/(3*(-x)^(3/10)) 本帖最后由 winner245 于 2014-7-14 23:09 编辑
刚刚对这个帖子里的公式用 latex 命令 （）重新编辑了一下，果然显示效果好多了，赞一个！
不过，对于发贴还是有些不尽如人意的地方，借此机会吐槽两句：
1. 图片显示限制：截图自动被限制在屏幕宽度的一半，这使得截图里的信息显得太小，而不易于浏览。当然点击图片后，图片会放大展开，但是如果图片多了，每一个都要单独点击，则很麻烦，不利于整体浏览。
2. 发帖验证限制：目前论坛对于主题帖的发表需要获取手机验证码才能完成。尽管这一限制的初衷可能是防止广告、灌水等滥发帖行为，但这一限制却可能对真正愿意在论坛分享东西的人造成很大的不便。如果发帖者在国内，这或许也算不上什么大问题，毕竟人手一部手机是很正常的事。可是，论坛应该也是面向国际的吧，海外那么多matlab爱好者也会上论坛，他们如果也要从一个注册地在国内的手机里获取验证码，那是何等的不便？这几次我发帖都需要向国内的亲朋好友求助，借助他们的手机才能完成发帖，真的觉得很麻烦。真心希望论坛能酌情考虑一下这一情况，使海外用户也能方便的在论坛上和大家交流
希望我的两句吐槽，不要引起什么误会，我只是一个单纯的matlab爱好者，单纯的和大家交流matlab心得
在将本文提到 2014a 符号积分 bug 反馈给 Mathworks 后，他们已经在最新的 2014b 版圆满解决了问题。测试如下：
>> syms x
>> double(int(x^2*exp(-x^2)*besseli(0,x),x,0,inf))
ans =
0.6380
>> syms a b positive
>> int(x^2*exp(-a*x^2)*besseli(0,b*x),x,0,inf)
ans =
(pi^(1/2)*exp(b^2/(8*a))*whittakerM(-1, 0, b^2/(4*a)))/(2*a*b)
其结果已经与 Maple 和 Mathematica 吻合 自从弃用了maple之后，matlab的符号运算就不敢再恭维了。现在已经安装maple。
页:
[1]