自适应辛普森法

我发现我纯属闲着没事干。

这玩意拿来求定积分。

辛普森公式

对于二次函数 \(f(x)=ax^2+bx+c\)，有：

\[\int_l^rf(x)\text dx=\frac{(r-l)(f(l)+f(r)+4f(\frac{l+r}2))}6 \]

证明大力拆。dirty work。

普通辛普森法

给个自然数 \(n\)，把区间 \([l,r]\) 分成 \(2n\) 个等长区间，每个区间用辛普森公式拟合积分值，然后加起来。

oi-wiki 说这玩意误差是

\[-\frac 1{90}(\frac{r-l}2)^5f^{(4)}(\xi) \]

其中 \(xi\) 是 \([l,r]\) 中的某个值。当然我们不用他。

自适应辛普森法

这玩意可以平衡时间和精度的限制。我们发现问题在怎么分段。

有一个简单明了的方法：如果这一段的图象和二次函数足够相似就停。

怎么判断？每次当前段带入公式求积分，然后从中点砍两半分别求积分。如果差值在精度误差内就停。

当然一般我们还要有一个最小迭代次数，迭代次数到了才停。

代码

double simpson(double l,double r){
    double mid=(l+r)/2;
    return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans,int dep){
    double mid=(l+r)/2;
    double fl=simpson(l,mid),fr=simpson(mid,r);
    if(abs(fl+fr-ans)<=15*eps&&dep<0)return fl+fr+(fl+fr-ans)/15;
    return asr(l,mid,eps/2,fl,dep-1)+asr(mid,r,eps/2,fr,dep-1);
}
double calc(double l,double r,double eps){
    return asr(l,r,eps,simpson(l,r),12);
}

模板 1：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
double a,b,c,d,L,R;
double f(double x){
    return (c*x+d)/(a*x+b);
}
double simpson(double l,double r){
    double mid=(l+r)/2;
    return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans){
    double mid=(l+r)/2;
    double fl=simpson(l,mid),fr=simpson(mid,r);
    if(abs(fl+fr-ans)<=15*eps)return fl+fr+(fl+fr-ans)/15;
    return asr(l,mid,eps/2,fl)+asr(mid,r,eps/2,fr);
}
double calc(double l,double r,double eps){
    return asr(l,r,eps,simpson(l,r));
}
int main(){
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
    double ans=calc(L,R,1e-6);
    printf("%.6lf\n",ans);
    return 0;
}

模板 2：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
double a;
double f(double x){
    return pow(x,a/x-x);
}
double simpson(double l,double r){
    double mid=(l+r)/2;
    return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans,int dep){
    double mid=(l+r)/2;
    double fl=simpson(l,mid),fr=simpson(mid,r);
    if(abs(fl+fr-ans)<=15*eps&&dep<0)return fl+fr+(fl+fr-ans)/15;
    return asr(l,mid,eps/2,fl,dep-1)+asr(mid,r,eps/2,fr,dep-1);
}
double calc(double l,double r,double eps){
    return asr(l,r,eps,simpson(l,r),12);
}
int main(){
    scanf("%lf",&a);
    if(a<0){
        puts("orz");return 0;
    }
    double ans=calc(1e-8,20,1e-8);
    printf("%.5lf\n",ans);
    return 0;
}