自适应辛普森法

Walking_Dead

2019-12-31 14:08:12

Personal

# 从入门到放弃 ## 什么是自适应辛普森法 $\ \ \ \ \ \ $ 震惊!!百度居然搜不到 ## 自适应辛普森法是什么 $\ \ \ \ \ \ $ 虽然这两个问题似乎是一样的,实际上是不同的两个问题 $\ \ \ \ \ \ $ 我们考虑需要求对于一个函数的 $$\int^R_L f(x),dx$$ $\ \ \ \ \ \ $ 对于这种问题,如果是一个具体的函数,我们当然可以用定积分来求,但是,我们能否有一个更加通用的算法呢? $\ \ \ \ \ \ $ 考虑一下我们的泰勒展开,我们一直用一个函数去逼近我们需要求到的函数,我们这里是不是可以用类似的方法呢? $\ \ \ \ \ \ $ 我们考虑对于一个区间,我们用一个二次函数去逼近它,我们假定我们的二次函数就叫做$g(x)=A\times x^2+B\times x+C$ $\ \ \ \ \ \ $ 我们知道对于$g(x)$的原函数为 $$G(x)=\frac{A}{3}\times x^3+\frac{B}{2}\times x^2+C\times x$$ $\ \ \ \ \ \ $ 然后,对于 $$\int_{a}^{b}f(x),dx\approx\int_{a}^{b}g(x),dx=G(b)-G(a)$$ $$=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(b-a)$$ $$=\frac{b-a}{6}(2A(a^2+ab+b^2)+3B(b+a)+6C)$$ $$=\frac{b-a}{6}(Aa^2+Ba+c+Ab^2+Bb+c+4A(\frac{a+b}{2})^2)+4B(\frac{a+b}{2})+4C)$$ $$=\frac{b-a}{6}(f(a)+f(b)+4f(\frac{a+b}{2}))$$ $\ \ \ \ \ \ $ 所以我们的就得到了: $$\int_{a}^{b}f(x),dx\approx\frac{b-a}{6}(f(a)+f(b)+4f(mid))$$ $\ \ \ \ \ \ $ 其中,$mid=\frac{a+b}{2}$ $\ \ \ \ \ \ $ 也就是说,对于一段区间我们有了一个大致的值,于是,我们就可以对于我们自己分出来的一段一段的区间求这一段区间的值,最后累加就好了 $\ \ \ \ \ \ $ 但是,我们如何分出区间呢?我们假设$simpson(l,r)=\frac{b-a}{6}(f(a)+f(b)+4f(mid))$,我们可以在使用分治法,如果$simpson(l,r)$与$simpson(l,mid)+simpson(mid,r)$大致相等,我们就可以不分了 ## 自适应辛普森法的模板 ```cpp double f (double x){return 为所求的函数;} double S (double l,double r) { double mid = (l + r) / 2; return (f (l) + 4 * f (mid) + f (r)) * (r - l) / 6; } double Abs (double x){return x > eps ? x : -x;} double Simpson (double L,double R) { double mid = (L + R) / 2; double Get1 = S (L,mid) + S (mid,R),Get2 = S (L,R); if (Abs (Get1 - Get2) <= 15 * eps) return Get1 + (Get1 - Get2) / 15; else return Simpson (L,mid) + Simpson (mid,R); } ``` ## 例题讲解 ### 【模板】自适应辛普森法1 [题目传送门](https://www.luogu.com.cn/problem/P4525) ### 思路 $\ \ \ \ \ \ $ 这道题可以用很多方法来做,比如[定积分(我的题解)](https://www.luogu.com.cn/blog/tanrui-2960967961/solution-p4525) $\ \ \ \ \ \ $ 但实际上就是一个模板 ### 代码 ```cpp #include <cmath> #include <cstdio> #include <cstring> #include <iostream> using namespace std; #define Int register int #define eps 1e-9 double a,b,c,d; double f (double x){return (c * x + d) / (a * x + b);} double S (double l,double r) { double mid = (l + r) / 2; return (f (l) + 4 * f (mid) + f (r)) * (r - l) / 6; } double Abs (double x){return x > eps ? x : -x;} double Simpson (double L,double R) { double mid = (L + R) / 2; double Get1 = S (L,mid) + S (mid,R),Get2 = S (L,R); if (Abs (Get1 - Get2) <= 15 * eps) return Get1 + (Get1 - Get2) / 15; else return Simpson (L,mid) + Simpson (mid,R); } signed main() { cin >> a >> b >> c >> d; double L,R; cin >> L >> R; printf ("%.6f\n",Simpson (L,R)); return 0; } ``` ### 【模板】自适应辛普森法2 [题目传送门](https://www.luogu.com.cn/problem/P4526) ### 思路 $\ \ \ \ \ \ $ 这道题看上去很复杂,因为里面有一个$\infty$,但是实际上,如果它是收敛的话,大概$x$在$12$左右就已经趋于$0$了,对答案不会再有影响,所以直接用自适应辛普森法高就好了 $\ \ \ \ \ \ $ 唯一的问题就是如何判断这个函数是否收敛了,我们可以看出,在$a<0$的时候这个函数显然是不收敛的,其他情况我们可以理性猜测都是收敛的 $\ \ \ \ \ \ $ 就...解决了? ### 代码 ```cpp #include <cmath> #include <cstdio> #include <cstring> #include <iostream> using namespace std; #define Int register int #define eps 1e-8 double a; double Abs (double x){return x > eps ? x : -x;} double f (double x){return pow (x,a / x - x);} double F (double l,double r) { double mid = (l + r) / 2; return (r - l) * (f(l) + 4 * f (mid) + f (r)) / 6; } double Simpson (double l,double r) { double mid = (l + r) / 2; double Get1 = F (l,mid) + F (mid,r),Get2 = F (l,r); if (Abs (Get1 - Get2) <= 15 * eps) return Get1 + (Get1 - Get2) / 15; else return Simpson (l,mid) + Simpson (mid,r); } signed main() { cin >> a; if (a < 0) puts ("orz"); else printf ("%.5f\n",Simpson (eps,11)); return 0; } ```