辛普森法

前置知识

定积分（少量），分治。
定积分就是求函数 \(f(x)\) 在区间 \([a,b]\) 中的图像包围的面积。（有正负，\(x\) 轴上为正，否则为负)。

应用范围

给你一个定积分求它的值。

思路

考虑小学的时候如何求一些不规则的图形的面积，我们可以划分成规则的图形求出面积在加起来。

我们常用的方法比如选两个点 \((x,f(x)),(y,f(y))\) 求之间的梯形面积，只要 \(x,y\) 之间隔得足够近，答案就是非常接近的。但是事实上时间通常允许分的足够小，所以误差较大。

我们考虑可以用曲线来拟合（误差可能可以小一点），所以我们思考用二次函数图像来逼近。

\[\begin{aligned} \int_{a}^{b} f(x)&\sim \int_{a}^{b}Ax^2+Bx+C\\ &={A\over3}b^3+{B\over2}b^2+Cb-{A\over3}a^3+{B\over2}a^2+Ca 微积分基本定理\\ &={A\over3}(b^3-a^3)+{B\over2}(b^2-a^2)+C(b-a)\\ &...一系列配凑\\ &={(b-a)\over 6}(f(a)+f(b)+4f({a+b\over 2})) \end{aligned} \]

详见
简单记忆，开头，结尾和 \(4\) 倍中间的函数值之和，除以系数和 \(6\) 乘上长度。

改进

我们需要平衡时间和精度，我们可以让代码自己根据精度要求来调整。
具体地，我们只需要算出 \([l,r],[l,mid],[mid,r]\) 区间的值判断大区间与小区间和之间的差值，如果已经达到精度要求就停止，否则递归进入 \([l,mid],[mid,r]\) 。正确性显然。

模板

double f(double x){return ...;}//将x代入定积分中的表达式
double si(double l,double r)
{
	double mid=(l+r)/2;
	return (r-l)*(f(l)+f(r)+4*f(mid))/6;//辛普森的式子
}
double getans(double l,double r)
{
	double mid=(l+r)/2;
	double zl=si(l,mid),zr=si(mid,r),z=si(l,r);
	if(fabs(zl+zr-z)<=eps)return zl+zr;//如果精度达到要求
	return getans(l,mid)+getans(mid,r);
}

根据思路来说是上面这样写，但是我们发现网上的模板有一个神秘的 \(15\) 。(貌似这样写会更加精准)

double getans(double l,double r,double eps)
{
	double mid=(l+r)/2;
	double zl=si(l,mid),zr=si(mid,r),z=si(l,r);
	if(fabs(zl+zr-z)<=15*eps)return zl+zr+(zl+zr-z)/15;
	return getans(l,mid,eps/2)+getans(mid,r,eps/2);
}

例题

模板1

将模板里的 \(f\) 补上即可。

点击查看代码

#include <bits/stdc++.h>
using namespace std;
double a,b,c,d,L,R,eps=1e-9;
double f(double x){return (c*x+d)/(a*x+b);}
double si(double l,double r)
{
	double mid=(l+r)/2;
	return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
double getans(double l,double r)
{
	double mid=(l+r)/2;
	double zl=si(l,mid),zr=si(mid,r),z=si(l,r);
	if(fabs(zl+zr-z)<=eps)return zl+zr;
	return getans(l,mid)+getans(mid,r);
}
int main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	cin>>a>>b>>c>>d>>L>>R;
	cout<<fixed<<setprecision(6)<<getans(L,R)<<'\n'; 
	return 0;
} 

模板2

难度上升。
首先没有保证收敛，可以通过一些方法证明 \(a<0\) 时是不收敛的，特判掉即可。(数学太差了，可以自行查阅资料证明)
然后区间范围有一个无限大，但是我们发现随着 \(x\) 的增大，\(f(x)\) 在急速的缩小，我们可以简单算一下在 \(a\) 取最大值时，\(x\) 在 \(10\) 附近就小于精度了。为了答案更精准我们可以区间再多算一点，算 \((0,15]\) 的答案（\(x=0\) 的函数无意义）。
然后就和模板一一样了。

点击查看代码

#include <bits/stdc++.h>
using namespace std;
double a,eps=1e-10;
double f(double x){return pow(x,(a/x)-x);}
double si(double l,double r)
{
	double mid=(l+r)/2;
	return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
double getans(double l,double r)
{
	double mid=(l+r)/2;
	double zl=si(l,mid),zr=si(mid,r),z=si(l,r);
	if(fabs(zl+zr-z)<=eps)return zl+zr;
	return getans(l,mid)+getans(mid,r);
}
int main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	cin>>a;
	if(a<0)cout<<"orz"<<'\n'; 
	else cout<<fixed<<setprecision(5)<<getans(eps,20.0)<<'\n'; 
	return 0;
}