matlab 做积分函数曲线

Question

函数如图：

想画出此函数曲线 r为横轴

自己写了一个 不对 求改正
x=0:0.01:15;
mu=0;b=1;s=1;
b0=b.^2;d0=s.^2;
f0=@(z)(1./z).*exp(-((log(z) - mu).^2./(2*d0))-(((x.^2)+(z.^2))./ (2*b0)));f1=quad(f0,0,inf);
f=(f1.*x)./(b0.*sqrt(2*pi*d0));
plot(x,f);

Answer 1

问题主要出现在f0这个函数。
 
在s比较小而x比较大的时候，指数函数的值在数值意义上为0（小于realmin=2.2251e-308），而Bessel函数的值为无穷大（大于realmax=1.7977e+308），导致出现0*Inf的情况，结果为非数NaN。被积函数的值一旦出现NaN，数值积分就会失败，返回的结果也是NaN，最后得到f是一串NaN，所以画不出曲线来。
 
就数值积分函数quad或integral而言，其内部包含了对边界点出现NaN值的处理，例如quad采取的做法是把边界点进行微调eps*(b-a)，以试图避免NaN。integral函数也有类似的处理。这也是类似下面的积分
1quad(@(x)sin(x)./x,0,1)
能够计算的原因（否则在x=0处会出现NaN）。
 
但数值积分函数对边界点所做的微调对于本题的被积函数是不够的。对此我花了不少时间对此进行分析，但现在有点太晚了，先不写了，把结果先给你吧。
 
可以把积分下限0改成0.1，这样会造成少量的误差，但影响不大：
123456789101112b=1; sigma=1; z=1; mu=0;b0=b.^2;d0=sigma.^2;k=(z.^2)./(2*b0);x=0:0.01:15;f0=@(s,x)(1./(s.^3)).*exp(-((log(s) - mu).^2./(2*d0))-(((k +1)*(x.^2))./(s.^2))).*besseli(0,(2*x*sqrt(k*(k+1)))./s); f1=arrayfun(@(x)integral(@(s)f0(s,x),0.1,inf),x);f=(2*(k+1)*exp(-k)*x.*f1)./(sqrt(2*pi*d0));plot(x,f,'r');
进一步的分析我明天贴到你对本问题的另外两处提问（编号1796169735918091867、1638082848257894860）。

Answer 2

1、quad的积分表达式中

2、quad的积分限不能为无穷大，换用integral函数（2012a以上版本）。

 

把

f0=@(z)(1./z).*exp(-((log(z) - mu).^2./(2*d0))-(((x.^2)+(z.^2))./ (2*b0)));
f1=quad(f0,0,inf);

改成

f0=@(z,x)(1./z).*exp(-((log(z) - mu).^2./(2*d0))-(((x.^2)+(z.^2))./ (2*b0)));
f1=arrayfun(@(x)integral(@(z)f0(z,x),0,inf),x);

本回答被提问者采纳