过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载,群号:927550334
MATLAB数值计算中,积分和微分是很重要的组成部分,积分可以梯形法做计算,也可以使用Monte Carlo方法做积分计算。过冷水在具有使用Monte Carlo方法做计算的时候取犯了一点小毛病,特此和大家分享,希望对大家精进Monte Carlo方法有帮助。
y=[0 0 0 0 0 0 0 0 1.71 3.24 4.59 5.76 6.75 7.56 8.19 8.64 8.91 9 8.91 8.64 8.19 7.56 6.75 5.76 4.59 3.24 1.71 0]'; x =[ -8 -7 -6 -5 -4.5 -4.3 -4.1 -4 -3.7 -3.4 -3.1 -2.8 -2.5 -2.2 -1.9 -1.6 -1.3 -1 -0.7 -0.4 -0.1 0.2 0.5 0.8 1.1 1.4 1.7 2 ]'; f= fit( x, y,'smoothingspline' ); figure1 = figure; axes1 = axes('Parent',figure1,'Position',[0.13 0.11 0.8 0.8]); hold(axes1,'on'); plot(x,y,'LineWidth',3,'Color',[1 0 0]); ylabel('$f(x)$','FontSize',24,'Interpreter','latex'); xlabel('$x$','FontSize',24,'Interpreter','latex'); xlim(axes1,[-10 6]); ylim(axes1,[-1 14]); hold(axes1,'off'); set(axes1,'FontSize',16); annotation(figure1,'arrow',[0.13 0.97],[0.11 0.11],'LineWidth',2); annotation(figure1,'arrow',[0.58 0.58],[0.11 0.95],'LineWidth',2);
该函数是一个简单函数用符号法、数值积分都可以获得精确解,实际案例函数很复杂,过冷水在此进行了模型简化,是为了说明问题,不要觉得过冷水在“杀鸡牛刀”。
需要用大数定理求积分
用大数定理求如图所示函数的具体表达式则为(N=1000)
1=linspace(-8,2,1000); y1=f(x1); I1=10*sum(y1)./1000 I1 = 35.9138
注意我们函数图像其实在(-8~-4)区间为0所以也可以等于:
x2=linspace(-4,2,1000); y2=f(x2); I2=6*sum(y2)./1000 I2 = 35.9806
从理论上分析没有问题,实际操作 it is error! 。这是过冷水犯的第一个错误。问题是时间需求是要从f(x)不为零的地方开始积分,过冷水为了省事做整个区间的积分,就是抱着对零积分不影响结果。
实际问题要稍微做做一下变形,要求分段积分。
当然如果你微积分还会一点点自然而然知道:
这个时候过冷水就沙雕了,用大数定理的或分别在(-8,-1),(-1,2),取了N(1000)个点,所有就有2N个点,所以应用大数定理就是:
x3=linspace(-8,-1,1000); x4=linspace(-1,2,1000); y3=f(x3);y4=f(x4); I3=10*sum([y3;y4])/(2000) I3 = 42.8187
再一次 it is error!到现在为止过冷水都没想明白为什么是错的,有对该问题感兴趣的可以留言讨论。“快就是慢,省事就是费事”简单问题还得从原理入手才可靠。
x3=linspace(-8,-1,1000); x4=linspace(-1,2,1000); y3=f(x3);y4=f(x4); I4=7*sum(y3)/1000 3*sum(y4)/1000 I4 = 35.9588
经过简单对比自然而然会发现
在早期过冷水做计算的时候一直使用的是I3,直到出现非常明显的错误的时候,一步一步做核查才知道是最基础的公式用错了,之前怀疑是自己代码什么地方敲错了。花了一晚上去检查错误最终排查完后非常有开心なるほど!
再回过头来讲讲为什么I1(35.9138)≠I2(35.9806)由于(-8~-4)函数值为零,在该区间取点分占一部分权重,(-8~-4)取点多了 (-4~2)取点自然就少一点 平均值就小一点。通常小数点后两位的误差我们是可以接受的。好巧不不巧过冷水在自己的问题计算误差放大的不合适,该错误可以用重要性抽样来解决。理论没问题,实际操作的时候就有问题了,还是不要省事严格按要求来。
这两处错误让我对积分和大数定理有了更深刻的了解,在简单的问题都要“勤实践,多思考”。总能从中学到新知识。本期的过冷水的分享就这么多,希望都能对大家有帮助,读者在实践中有遇到相似的问题,都可以留言分享,大家共同学习进步。
过冷水发表于 仿真秀 平台原创文章,未经授权禁止私自转载,如需转载请需要和作者沟通表明授权声明,未授权文章皆视为侵权行为,必将追责。如果您希望加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载均可加群:927550334。
精品回顾
过冷水和你分享 matlab读取存储各种文件的方法 文末有独家金曲分享