matlab画积分函数曲线
b=1; h=log(10)/20;s=h*b;z=1; u=0;b0=b.^2;d0=s.^2;k=(z.^2)./(2*b0);x=0:0.01:15;
f0=@(v,x)(1./(v.^3)).*exp(-(((log(v) - u).^2)./(2*d0))-(((k +1)*(x.^2))./(v.^2))).*besseli(0,(2*x*sqrt(k*(k+1)))./v);f1=integral(@(v)f0(v,x),0,inf); f2=arrayfun(@(x)f1,x);
f=(2*(k+1)*exp(-k)*x.*f2)./(sqrt(2*pi*d0));plot(x,y,'r'); 展开
这个问题和另一个问题(编号2051722037141864067)基本相同,但与那个问题相比,又多了一处错误:
f1=integral(@(v)f0(v,x),0,inf);
f2=arrayfun(@(x)f1,x);
和
f2=arrayfun(@(x)integral(@(v)f0(v,x),0,inf),x);
的作用是不一样的:
对于后一种写法,在integral的表达式中,x是作为匿名函数的形式参数,其值不定,在调用arrayfun的时候才被赋予了具体的值,即向量x=0:0.01:15,也就是说,在这个表达式里出现的x有两种不同的含义,其实是两个不同的量使用了相同的符号,这一点初学者有点容易被搞糊涂。
而对于前一种写法,f1不是匿名函数,而是直接把x=0:0.01:15代入求得函数值,但由于计算积分的时候,自变量v是按照向量处理的,而此时x也是向量,所以会导致出错,无法计算。
在另一个问题(编号2051722037141864067)的回答里面,我分析了造成无法绘图的原因是由于被积分函数的值出现了NaN,并建议把积分的下限改为0.1。这样并不会造成太大影响,对此,我们可以通过下面的方法说明:取积分下限分别为0.1、0.01、0.001、0.0001,分别计算并绘图(匿名函数f0之前的代码基本不变,但按照你上一次提问的值,取s=1):
c = 'rgbm';
for i=1:length(c)
f2=arrayfun(@(x)integral(@(v)f0(v,x),10^-i,inf),x);
f=(2*(k+1)*exp(-k)*x.*f2)./(sqrt(2*pi*d0));
plot(x,f,c(i));
hold on
end
legend(arrayfun(@(i){sprintf('积分下限 = %g',i)},10.^-(1:4)),0)
可以看到,积分下限取0.1对于x=0~15都能够计算,取0.01则只能计算大约0~4的范围(x更大的时候会得到NaN,那部分的曲线就画不出来了),如果进行局部放大:
可以看到,如果取积分下限为0.001和0.0001,所给积分表达式能够正确计算的范围分别只有0.4和0.04左右。另一方面,我们注意到,在可计算的范围内,积分下限取0.01、0.001、0.0001的曲线是相同的,这也就意味着,积分下限取0.01就足够精确了。另外,从图中也反映出,积分下限如果取0.1,会在大约x=0~0.2的范围内存在误差。
如果按照你现在的s取值,画出的图如下,与s=1相比,取相同的积分下限,可计算的范围总体上相同,但即使取0.1作为积分下限,也几乎没有误差。
另外,还有一些小问题:
1、计算出来的结果是f,绘制曲线却用y。
2、括号并不是越多越好,过多的括号反而看起来更乱,例如
f0=@(v,x)(1./(v.^3)).*exp(-(((log(v) - u).^2)./(2*d0))-(((k +1)*(x.^2))./(v.^2))).*besseli(0,(2*x*sqrt(k*(k+1)))./v);
可以改成
f0=@(v,x) (1./(v.^3)) .* exp(-((log(v) - u).^2./(2*d0))-(k +1)*x.^2./v.^2) .* besseli(0,2*x*sqrt(k*(k+1))./v);
看起来更清晰。