自适应辛普森了解一下

Posted 2021-12-07 judge

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