A 了这道超短的紫题,来发表一下自己的一些想法...
简单介绍辛普森这玩意儿
不如先学学泰勒展开?
首先泰勒展开大家都听说过吧?【雾
没听说过?安利某知乎回答:苍老师教你如何更好地记忆泰勒展开
然后你就知道了,泰勒展开其实是对于某个函数在一个点不断去高阶求导,然后用求导得到的信息构造一个多项式,使得这个多项式在一定范围内几乎和原函数拟合(可以理解为接近重合的意思吧...)
类比到辛普森?
那么其实自适应辛普森也是类似的道理,只不过它是用了分治的方法去构造这个 假 拟合 多项式(其实就是二次函数,究其原因应该就是这玩意儿比较好积分并且存在弧度,容易拟合吧...)
自适应是个什么鬼?
至于为什么前面加了个自适应呢?因为我们考虑分治是要有终止条件的(像对于一个序列分治的话就是分治的区间长度为 1 时停止),但是我们这里是在实数域上拟合一个多项式啊,不存在什么规定的终止点...
于是我们考虑怎样去设置终止条件?那当然是给定一个误差范围(比如 1e-6),然后如果拟合函数和原函数大多数点的 y 值误差不超过这个范围就终止分治,直接拿当前的函数去积分就好了,当然具体怎么比较的先别管(你看到下面之后会发现根本不需要比较两个函数 2333 )
然后我们发现这个误差越大答案越不准确,越小跑的越慢,那么我们就要调整这个误差范围,然后这样的过程就有点像"自适应"了
拟合函数怎么构造?
我们先将原函数约等于成一个二次多项式:
\[f(x)≈Ax^2+Bx+C\]
然后题目要我们求的东西也转化一下:
\[ANS=\int_a^b f(x)dx\]
\[≈\int_a^b Ax^2+Bx+C\]
\[={A\over 3}(b^3-a^3) + {B\over 2} (b^2-a^2) +C(b-a)\]
\[={(b-a)\over 6} [ 2A(b^2+ab+a^2) + 3B(b+a)+6C]\]
\[={(b-a)\over 6} [ (Aa^2 + Ba+C) + (Ab^2+Bb+C)+4(~A~({a+b\over 2})^2 + B ({a+b\over2})+C)]\]
\[≈{(b-a)\over 6} [~ f(a) + f(b) + 4~f({a+b\over2})~]\]
恩?你问我推导哪里来的? 大佬%%%
我们发现最后约回去了,A B C 都不见了,而且这时候式子里面已经没有积分了...
所以怎么构造 A、B、C 还是没讲?
没错,我们发现上面其实根本不需要用到这个拟合二次函数的具体系数 A、B、C ,只需要将 f 的式子带入计算求值就好了
但是这样的话我们发现之前说的终止条件不见了,因为我们的拟合函数已经不需要了(而且找出来也很麻烦...)
其实我们只需要比较当前区间 \([a,b]\) 不分治时候的答案 \(ans\) 和 分治下的答案 \(ansL+ansR\) 的误差是否超过了我们设置的精度范围就好了
至于正确性?(我怎么知道)
因为我们知道这个分治肯定是层数越多越精确的,所以分治的结果精确度肯定高于当前 ans 的精确度,那么我们卡卡精度就能让答案在误差允许范围内了
FAQ: mmp 原来真的不用比较原函数和拟合函数...
如果要比较的话又能怎么比较呢?反正我是不会的...
code
代码超级短!相应的我们可以知道这里我们只需要将 f 函数改一改就能解决其他函数的积分问题了(一定精度下)
//by Judge
#include<bits/stdc++.h>
#define db double
using namespace std;
db a,b,c,d,l,r;
inline db f(db x){return (c*x+d)/(a*x+b);}
inline db simpson(db l,db r){ db mid=(l+r)/2;
return (f(l)+f(r)+4*f(mid))*(r-l)/6;
}
db asr(db l,db r,db eps,db ans){
db mid=(l+r)/2,L=simpson(l,mid),R=simpson(mid,r);
if(fabs(L+R-ans)<=eps) return L+R+(L+R-ans);
return asr(l,mid,eps/2,L)+asr(mid,r,eps/2,R);
}
inline db asr(db l,db r,db eps){
return asr(l,r,eps,simpson(l,r));
}
int main(){ scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
return !printf("%.6lf\n",asr(l,r,1e-8));
}
原文地址:https://www.cnblogs.com/Judge/p/10927547.html