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);

来源:学生作业帮助网 编辑:作业帮 时间:2024/05/06 16:45:54
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);

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);
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;

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);

问题主要出现在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).