本帖最后由 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(数学符号记作: ,MATLAB 中对应的函数是 besseli),介绍如何有效的克服积分结果为 NaN 的问题。本帖讲述的方法虽然只是针对这一种特殊函数,但解决问题的思路适用于这一类数值计算问题,希望能对有兴趣的人有帮助!
最后,在分享这个案例的过程中,我还会顺带指出 MATLAB/MuPad 符号积分系统的一个计算错误,该错误已在 MATLAB R2014a 下验证过,应该算是一个 bug。
数值积分案例
给出我要计算的积分:
下面,我们来看看直接用 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 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数值计算
>> f=@(x)x^2.*exp(-x^2).*besseli(0,x) >> f(1000)
ans = NaN
数值计算得到结果为NaN。这正是因为单独的两项besseli(0,x) 和exp(-x^2)分别计算时得到了Inf和0,如下所示: ans = 0 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。积分式如下:
(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 了把。
|