matlab画积分函数曲线

说明

这个问题和另外两个问题（编号2051722037141864067、1638082848257894860）基本上是重复的，我已经在那两个问题做了回答，主要原因是匿名函数f0在x比较大、s比较小的情况下会出现非数NaN，导致计算失败。更具体的分析与解决方法这里不再重复，感兴趣的请自行查阅（因度娘经常抽风，就不贴链接了，把编号的数字复制替换本问题的地址中question后面的那一串数字即可）。

这里再对两个问题做进一步探讨：一是对出现NaN的原因做更深入分析，二是把积分下限换成0.1的误差分析。

1、结果中出现NaN的原因

之桐祥前分析过（参见问题2051722037141864067），之所以不能画图，归根到底是由于f0在某些条件下计算结果出现NaN引起的，而NaN又是由于指数项为0、Bessel函数为无穷大导致的。

■ 对于指数函数exp(-x)，在什么条件下结果为0？

负指数函数是x的减函数，从数学的角度来说，函数值会逐渐衰减趋近于0，但只要x是有限值，函数值就不真正为0。但从数值在计算机内的表示来说，双精度浮点数只有8个字节，其表示精度与范围都是有限的，可以判断，x为某个有限值的时候，函数值就会小于最小的正浮点数，也就是数值意义上0。

最小的正浮点数可以用realmin获得，但请注意，这个是所谓规格化（normalized）浮点数，而不是最小的浮点数，最小的浮点数是eps(realmin)（或eps*realmin）：

>> realmin
ans =
    2.225073858507201e-308
>> eps(realmin)
ans =
    4.940656458412465e-324

所谓normalized，是指大于该浮点数的运算能够保证精度，一旦小于realmin，被称为IEEE "denormal"，不能再保证运算精度。例如：

>> eps(realmin)/2
ans =
     0
>> eps(realmin)/1.99
ans =
    4.940656458412465e-324

我们看到，这个最小的浮点数除以1.99仍然等于其自身，除以2则等于0。事实上，这个数的浮点数表达只有最后一个bit是1，其它63bit都是0，一旦除以2或更大的数，就会得到全0的八个字节，也就是0。

了解了最小的浮点数，也就可以知道使得exp(-x)数值上达到0的x值了：

>> x=-log(eps(realmin))+log(2)
x =
  745.1332
>> exp(-x)
ans =
  4.9407e-324
>> exp(-745.1333)
ans =
     0

也就是说，这个数稍大于745。

■ Bessel函数什么条件下为无穷大？

这个有点遗憾，由于Bessel函数不像exp那样有逆函数可用，使得besseli(0,x)为无穷大的x我只能通过试探大致确定在700-701之间：

>> besseli(0,700)
ans =
  1.5296e+302
>> besseli(0,701)
ans =
   Inf

当然，可以通过进一步的试探确定更多的有效数字：

>> besseli(0,700.9217936944459)
ans =
  3.8426e+302 
>> besseli(0,700.921793694446)
ans =
   Inf

■ 函数避免出现NaN的条件

指数项包括两部分，仅以其中一项为例（未包含负号）：

ezplot(@(v)(log(v) - u).^2./(2*d0),[0 0.1])

可以看到，在v比较小时，仅此一项就会让负指数函数的值为0。另外一项对应的指数也是负数，对应的昌正函数值小于1。可以通过下面的方法大致看到

ezmesh(@(v,x) -((log(v) - u).^2./(2*d0))-(k +1)*x.^2./v.^2,[0 0.1 0 15])
zlim([-750 0])
view(0,90)

空白区域即意味着指数小于-750，也就是函数值为0。可见，指数项除了在小部分区域外，耐轮悔大多数条件下的函数值为0，这样，为了避免出现0*Inf，重点在于防止Bessel函数出现无穷大的值。而即使指数部分不为0，一旦Bessel函数值为Inf，两项相乘的值为Inf，计算结果同样没有意义。

看一下Bessel函数的变量：

ezplot(@(v,x) 2*x*sqrt(k*(k+1))./v - 700,[0 0.1 0 50]) 
axis auto

图中曲线的含义是，当v取某个值的时候，x只有小于特定值，Bessel函数才为有限值。这个值大约是对应v=0.01,x=4.04；v=0.1，x=40.4，也就是说，对于v=0.01，只能计算大约x<4.04范围的Bessel函数，v=0.1时，可计算范围大约是x<40.4。

这个结论和之前的分析吻合。另一方面，我们可以看到，如果只需要画x=0~15区间的积分函数，可以取积分下限为0.04。

2、把积分下限换成0.1的误差分析

按照问题1638082848257894860的分析，把积分下限进行微调成0.1，对于大部分的函数值没有影响。现在具体看看误差有多大。

由于只是对积分限进行微调，所以需要考虑的只是被积变量v在0~0.1区间f0函数的情况。这里按照sigma=1来分析（如果按照本题的s分析，几乎没有误差，这一点也可以在问题1638082848257894860里面的曲线看到）。

画出x取不同值的f0-v曲线：

ezplot(@(v)f0(v,0),[0 0.1])
hold on
ezplot(@(v)f0(v,0.01),[0 0.1])
ezplot(@(v)f0(v,0.05),[0 0.1])
ezplot(@(v)f0(v,0.1),[0 0.1])
axis auto

可以看到，x=0对应的曲线值是最大的（应该可以从理论上证明），但最大值也是有界的。对积分下限进行微调导致的误差不会超过这部分积分再乘一个相应的系数。

限于时间精力，这部分的分析未能进一步深入。写了也没几个人看，就先这样吧。

最后，对积分下限取0.01和0.1的误差进行比较：

x = 0:0.01:0.5;
df = arrayfun(@(x)integral(@(v)f0(v,x),0.01,inf),x) - arrayfun(@(x)integral(@(v)f0(v,x),0.1,inf),x);
f = (2*(k+1)*exp(-k)*x.*df)./(sqrt(2*pi*d0));
plot(x,f,'r');

这也和之前分析的吻合，即只在x比较小的区间（大约0.2）才有一定误差。而取0.001和0.01的误差更小：

结束语

花费好几个小时写的分析，很大程度上和解决楼主所提问题本身已经没有太大关系，只是为了探究使用MATLAB可能出现的误差或异常现象的深层原因，以便在以后的应用中加以注意。

在此，向楼主提个请求，能否告知这个积分函数的应用背景？花了这么多时间研究这个问题，虽然只是出于个人的爱好，但我把这些拿出来分享的时候，希望能够知道这究竟是哪个领域的问题，谢谢。