学习笔记: 自适应Simpson积分+例题

参考bloghttps://blog.csdn.net/lunch__/article/details/82288676
https://blog.csdn.net/xyz32768/article/details/81392369

这个东西对于三次以下的函数是正确的，但是对于三次以上的函数我们可以将其近似为低次函数套用Simpson公式，这个东西学名好像叫自适应Simpson积分。

昨天ACMACM模拟的时候遇到了一道SimpsonSimpson积分相关的题

完全不知道怎么求，我们组FishmenFishmen被BymBym嘲讽了很久

于是今天下午结合各种资料还是看了一下

这个东西不要觉得它看上去讲什么积分很高级 实际上认真推导也不是很难

首先SimpsonSimpson积分的式子是这个东西

积分号的下标代表区间左端点，上标代表区间右端点

这个值的意思就是函数在[a,b][a,b]区间与xx轴，x=ax=a， x=bx=b围成的面积

首先引入一些东西，SimpsonSimpson积分就是用一个二次曲线来无限逼近原曲线

来达到求得近似值的效果

首先引入一个东西 对于一个二次函数f(x)=ax2+bx+cf(x)=ax2+bx+c
求积分可得F(x)=∫x0f(x)dx=a3x3+b2x2+cx+DF(x)=∫0xf(x)dx=a3x3+b2x2+cx+D
在这里DD是一个常数 这个东西具体怎么证等我问了物理组的同学再来填坑

那么∫RLf(x)dx=F®−F(L)∫LRf(x)dx=F®−F(L)
∫RLf(x)dx=a3(R3−L3)+b2(R2−L2)+c(R−L)∫LRf(x)dx=a3(R3−L3)+b2(R2−L2)+c(R−L)
∫RLf(x)dx=(R−L)[a3(L2+R2+LR)+b2(L+R)+c]∫LRf(x)dx=(R−L)[a3(L2+R2+LR)+b2(L+R)+c]
∫RLf(x)dx=b−a6(2aL2+2aR2+2aLR+3bL+3bR+6c)∫LRf(x)dx=b−a6(2aL2+2aR2+2aLR+3bL+3bR+6c)
∫RLf(x)dx=b−a6{(aL2+bL+c)+(aR2+bR+c)+4[a(L+R2)2+b(L+R2)+c]}∫LRf(x)dx=b−a6{(aL2+bL+c)+(aR2+bR+c)+4[a(L+R2)2+b(L+R2)+c]}
∫RLf(x)dx=b−a6[f(L)+f®+4f(L+R2)]∫LRf(x)dx=b−a6[f(L)+f®+4f(L+R2)]
就这样证完了… 然后我们来考虑这个代码怎么实现

double Simpson(double a, double b) {
int mid = (a + b) / 2;
return (b - a) / 6 * (f(a) + f(b) + 4 * f(mid));
}
1
2
3
4
好了写完了
这样子直接计算我们发现误差非常大，毕竟原图像可能不能很好的用二次函数逼近

自适应Simpson积分
那怎么办呢，我们可以考虑模拟微分的过程，把区间分成无穷小份再把得到的答案积分就是最后的答案了

但是一般我们只要求一个近似值，所以我们就递归到满足正常误差即可

还有一个问题 每次递归下去都有精度误差，那么累加起来可能就和正确答案相差甚远了

这个时候我们可以在递归的时候提高精度要求，尽量让小的答案精确这样就可以了

double DFS(double L,double R,double eps){
    double ans=J(L,R);
    double mid=(L+R)/2;
    double ans1=J(L,mid),ans2=J(mid,R);
    if(fabs(ans1+ans2-ans)<eps*15) return ans1+ans2+(ans1+ans2-ans)/15;
    return DFS(L,mid,eps/2)+DFS(mid,R,eps/2);
}

在上式中eps一般取1e-6（具体看题题目的要求了，在下面的一道洛谷模板题中要求精度达到小数点后六位）

洛谷模板题
下面是AC的代码

//#include <cstdio>
//#include<cstring>
//#include<iostream>
//#include<algorithm>
//#include<set>
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,b) for(int i=a;i<=b;i++)
const int maxn = 0x3f3f3f3f;
const int N=1e6+10;
const double eps=1e-6;
double a,b,c,d,l,r;
double F(double x){
    return (c*x+d)/(a*x+b);
}
double J(double L,double R){
    double mid=(L+R)/2;
    return (R-L)/6*(F(R)+4*F(mid)+F(L));
}
double DFS(double L,double R,double eps){
    double ans=J(L,R);
    double mid=(L+R)/2;
    double ans1=J(L,mid),ans2=J(mid,R);
    if(fabs(ans1+ans2-ans)<eps*15) return ans1+ans2+(ans1+ans2-ans)/15;
    return DFS(L,mid,eps/2)+DFS(mid,R,eps/2);
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif // ONLINE_JUDGE
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6f\n",DFS(l,r,eps));
}

然后还有一道HDU1724的板子题
题意：题目就是求两条线所夹的椭圆的面积。这是个标准的椭圆，其中心在原点。我好像看到了题目连椭圆面积公式都给了，然而并没有任何用，直接上辛普森！
如果改变eps的大小，比如这道题eps=1e-6时46ms，eps=1e-4时15ms，看题目而定

//#include <cstdio>
//#include<cstring>
//#include<iostream>
//#include<algorithm>
//#include<set>
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,b) for(int i=a;i<=b;i++)
const int maxn = 0x3f3f3f3f;
const int N=1e6+10;
const double eps=1e-6;
double a,b;
double F(double x){
    return sqrt(b*b*(1-(x*x)/(a*a)));
}
double J(double L,double R){
    double mid=(L+R)/2;
    return (R-L)/6*(F(R)+F(L)+4*F(mid));
}

double DFS(double L,double R,double eps){
    double mid=(L+R)/2;
    double ans=J(L,R),ans1=J(L,mid),ans2=J(mid,R);
    if((ans1+ans2-ans)<=eps*15) return ans1+ans2+(ans1+ans2-ans)/15;
    return DFS(L,mid,eps/2.0)+DFS(mid,R,eps/2.0);
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif // ONLINE_JUDGE
    double l,r;
    int t;
    scanf("%d",&t);
    while(t--){
    scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
    printf("%.3f\n",2.0*DFS(l,r,eps));
    }
}