辛普森法

自适应辛普森法

引入：

对于计算一个非容易计算的多边形时，我们通常的想法就是在误差范围内将其划分成多个矩形，然后依次计算。

但是，如何划分这些矩形合适呢，这就是自适应辛普森法解决的问题。

方法：

简而言之，就是一种用二次函数来逼近被积函数。
把求原来函数的积分换成求二次函数的积分的一种近似求积分的方法

首先有公式：

\[\int_a^b f(x) dx \approx \frac{(b-a)}{6}(f(a)+f(b)+4f(\frac{a+b}{2})) \]

为什么叫自适应呢，其实就是以二分的方式把函数划分成一个个小的区间进行求和。

解决过程

具体就是分别算两边 \(simps(l,mid)\) 与 \(simps(mid,r)\)，比较两者的和与 \(simps(l,r)\) 的关系:

1. 如果误差足够小的话，就可以返回.
2. 否则继续二分计算。

例题：
P4525 【模板】自适应辛普森法1

将公式中的 \(a\) 和 \(b\) 分别转化为两个函数之间对应 \(a\) 位置和 \(b\) 位置的高度差，就可以直接带入公式计算。
代码：

#include<bits/stdc++.h>
using namespace std;
#define db double
const db eps=1e-8;
db A,B,C,D,L,R;
double hei(double x){return (C*x+D)/(A*x+B);}//计算高度之差
double simps(db a,db b,db fa,db fb){//计算三点辛普森公式的值，缩小矩形的面积
    db c=a+(b-a)/2,fc=hei(c);
    return (fa+4*fc+fb)*(b-a)/6.0;
}
/*
曲线下方划分成若干矩形，用矩形面积和估算曲线下方面积
容易近似的地方少划矩形，不容易近似多画矩形
*/
db asr(db a,db b,db fa,db fb,db eps){
    db c=a+(b-a)/2,fc=hei(c);
    db l=simps(a,c,fa,fc),r=simps(c,b,fc,fb),all=simps(a,b,fa,fb);
    if(fabs(l+r-all)<=15*eps) //精度充足，不需要分
        return l+r+(l+r-all)/15.0;
    return asr(a,c,fa,fc,eps/2)+asr(c,b,fc,fb,eps/2);//精度不充足
}
db solve(db a,db b){return asr(a,b,hei(a),hei(b),eps);}
int main()
{
    cin>>A>>B>>C>>D>>L>>R;
    printf("%.6lf\n",solve(L,R));
    system("pause");
    return 0;
}

P4526 【模板】自适应辛普森法2
根据积分的性质，发散就是 \(a<0\)

然后进行计算就行

#include<bits/stdc++.h>
using namespace std;
#define db double
const double eps=1e-9;
db a,l,r,mid;
inline db f(db x){return pow(x,a/x-x);}
inline db simps(double l,double r){
    return (r-l)*(f(l)+f(r)+4.0*f((l+r)/2.0))/6.0;
}
inline db solve(double l,double r,double ans){
    double mid=(l+r)/2.0,L=simps(l,mid),R=simps(mid,r);
    if(fabs(L+R-ans)<=eps) return  ans;
    return solve(l,mid,L)+solve(mid,r,R);
}
int main()
{
    cin>>a; if(a<0) cout<<"orz"<<endl;
    else printf("%.5lf\n",solve(eps,20,simps(1e-9,20)));
    system("pause");
    return 0;
}