Question

matlab 做积分函数曲线

想画这条曲线 编程如下：
b=1; s=1; z=1; u=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,inf),x);
f=(2*(k+1)*exp(-k)*x.*f1)./(sqrt(2*pi*d0));
plot(x,y,'r');

可是有很多warning,而且出不来曲线 不知道问题在哪里 谢谢帮忙
sigma=s;
mu=u;

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函数也有类似的处理。这也是类似下面的积分

quad(@(x)sin(x)./x,0,1)

能够计算的原因（否则在x=0处会出现NaN）。

 

但数值积分函数对边界点所做的微调对于本题的被积函数是不够的。对此我花了不少时间对此进行分析，但现在有点太晚了，先不写了，把结果先给你吧。

 

可以把积分下限0改成0.1，这样会造成少量的误差，但影响不大：

b=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

sigma为多少？mu为多少？

追问

Oh sorry sigma=s mu=u

追答

是一条直线吗？

1.m

b=1; z=1; sigma=1;mu=0;b0=b.^2;d0=sigma.^2;k=(z.^2)./(2*b0);
syms x
eq1=2*(k+1).*x/sqrt(2*pi*d0)*exp(-k);
eq2=1./sigma^3.*exp(-(k+1)*x.^2)/sigma^2-(log(sigma)-mu)^2/(2*d0).*besseli(mu,2*x*sqrt(k*(k+1))./sigma);
y=eq1*int(eq2,0,inf)

2.m

x=0:0.1:15;
y=2731571253070597/22577700327169044.*x*3^(1/2)*2^(1/2)*pi^(1/2);

plot(x,y,'-','LineWidth',2),xlabel('r'),ylabel('f(r)')