过冷水最近这段时间在做积分的学习工作,之前连续分享了好几期的 蒙特卡洛法应用;你所不知道的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)
x1=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)取点自然就少一点 平均值就小一点。通常小数点后两位的误差我们是可以接受的。好巧不不巧过冷水在自己的问题计算误差放大的不合适,该错误可以用重要性抽样来解决。理论没问题,实际操作的时候就有问题了,还是不要省事严格按要求来。
这两处错误让我对积分和大数定理有了更深刻的了解,在简单的问题都要“勤实践,多思考”。总能从中学到新知识。本期的过冷水的分享就这么多,希望都能对大家有帮助,读者在实践中有遇到相似的问题,都可以留言分享,大家共同学习进步。