自适应辛普森算法

自适应辛普森公式

 

名字很高大上，事实上是计算机与数学深度结合后诞生的一种算法。

求积分，我们知道一定存在一个函数，可以表示这段区间图像。

辛普森公式（请读者自行百度）：

基本思想是我们把三个点看做二次函数的一部分，用二次函数的拟合原本函数的积分。

对于一般的二次函数公式，我们可以推出：

 

采用极限的思路，我们可以不断加大精度。

事实上，我们并不需要其中微不足道的那部分，只需保证我们所需的精度即可

 辛普森公式函数：

int simpson(int l, int r) { //这里做了 #define double int
    // f(x) 为题目要求的函数
    int mid = (l + r) / 2;
    return (r - l) * (f(l) + f(r) + 4 * f(mid)) / 6;
}

 

那么什么是自适应呢？就是控制精度的过程。不断递归，直到误差处于可以忽略的范围。

假设我们要求 [l,r] 上的积分，则先用辛普森公式求出这一段用抛物线拟合的积分，再分成两段递归求积分并加和。如果这加和结果和拟合出来的结果相差很小，就可以直接跳出递归。

 

int adaptiveSimpson(int l, int r, int eps, int res, int dep) {
    double mid = (l + r) / 2, L = simpson(l, mid), R = simpson(mid, r);
    if (abs(L + R - res) <= 15 * eps && dep < 0)
        return L + R + (L + R - res) / 15;
    return adaptiveSimpson(l, mid, eps / 2, L, dep - 1) +
        adaptiveSimpson(mid, r, eps / 2, R, dep - 1);
}

int calc(int l, int r, int eps) {
    return adaptiveSimpson(l, r, eps, simpson(l, r), 12);
}

 

 eg. P4525 【模板】自适应辛普森法 1

 代码

#include<bits/stdc++.h>
#define int long double
using namespace std;
const int N=1e6+10;
int a,b,c,d,l,r;
int f(int x){
    return (c*x+d)/(a*x+b);
}
int simpson(int l,int r){
    int mid=(l+r)/2;
    return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
int adaptiveSimpson(int l,int r,int eps,int res,int dep){
    double mid=(l+r)/2,L=simpson(l,mid),R=simpson(mid,r);
    if(abs(L+R-res)<=15*eps&&dep<0)
        return L+R+(L+R-res)/15;
    return adaptiveSimpson(l,mid,eps/2,L,dep-1)+
        adaptiveSimpson(mid,r,eps/2,R,dep-1);
}
int calc(int l,int r,int eps){
    return adaptiveSimpson(l,r,eps,simpson(l,r),12);
}
signed main(){
    cin>>a>>b>>c>>d>>l>>r;
    printf("%.6Lf",calc(l,r,1e-6));
    return 0;
}