自适应辛普森算法

一个并不知道现在学了有什么用的算法，因为我不会计算几何，甚至是对计算几何一窍不通。

自适应辛普森算法（ASR）可以用来求 \(\int_a^bf(x)\,\mathrm dx\)，即求一个函数 \(f(x)\) 的定积分，其中 \(f(x)\) 是一个不太好直接积的函数。其大致思路大概就是不断用一个二次函数对原函数进行拟合，然后用二次函数的积分方式对原函数求近似积分，并通过左右子区间的辛普森值来判断误差，如果误差在精度范围内就直接返回拟合得到的二次函数的积分值。

自适应辛普森算法的适用前提是 \(f(x)\) 为平滑的函数，如果出现下图的情况那就凉凉了（不过一般出题人也不会恶意卡？）

Simpson 公式的推导

我们考虑对 \(f(x)\) 进行拟合，假设 \(f(x)\approx Ax^2+Bx+C\)

那么：

\[\begin{aligned} &\int_a^bf(x)\,\mathrm dx\\ \approx&\int_a^bAx^2+Bx+C\\ =&\dfrac{A}{3}(b^3-a^3)+\dfrac{B}{2}(b^2-a^2)+C(b-a)\\ =&\dfrac{A}{3}(b-a)(b^2+ab+a^2)+\dfrac{B}{2}(b-a)(b+a)+C(b-a)\\ =&\dfrac{b-a}{6}(2A(b^2+ab+a^2)+3B(a+b)+6C)\\ =&\dfrac{b-a}{6}(2Ab^2+2Aab+2Aa^2+3Ba+3Bb+6C)\\ =&\dfrac{b-a}{6}((Aa^2+Ba+C)+(Ab^2+Bb+C)+(Aa^2+2Aab+Ab^2+2Ba+2Bb+4C))\\ =&\dfrac{b-a}{6}((Aa^2+Ba+C)+(Ab^2+Bb+C)+(4A(\dfrac{a+b}{2})^2)+4B(\dfrac{a+b}{2})+4C)\\ =&\dfrac{b-a}{6}(f(a)+f(b)+4f(\dfrac{a+b}{2})) \end{aligned} \]

也就是说 \(\int_a^bf(x)\,\mathrm dx\) 可以用 \(\dfrac{b-a}{6}(f(a)+f(b)+4f(\dfrac{a+b}{2}))\) 来近似地求出。

---

自适应辛普森算法大概也就是这个思路，不过直接套公式显然精度误差太大，我们考虑这样的算法：我们要求 \([l,r]\) 内的积分值，我们先用辛普森公式算出 \([l,r]\) 的近似值 \(A\)，再设 \(\dfrac{l+r}{2}=mid\)，计算出 \([l,mid],[mid,r]\) 的近似值 \(L,R\)，如果 \(|L+R-A|\) 在误差范围内就返回 \(L+R+(L+R-A)\)，否则继续递归，这也就能得到近似值了。

复杂度 \(\mathcal O(\text{玄学})\)，正确性玄学（大雾）

因为我不会计算几何所以也没啥应用可刷，就刷刷模板题罢：

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

模板题，直接套用上面的算法即可。

const double EPS=1e-9;
double a,b,c,d,l,r;
double func(double x){return 1.0*(c*x+d)/(a*x+b);}
double simpson(double l,double r){
	double mid=(l+r)/2;
	return (func(l)+4*func(mid)+func(r))*(r-l)/6;
}
double inter(double l,double r,double A){
	double mid=(l+r)/2;
	double L=simpson(l,mid),R=simpson(mid,r);
	if(fabs(L+R-A)<=EPS) return L+R+(L+R-A);
	else return inter(l,mid,L)+inter(mid,r,R);
}
int main(){
	scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
	printf("%.6lf\n",inter(l,r,simpson(l,r)));
	return 0;
}

似乎也有纯数学的推导，赶紧补一下，正好强化下我微积分的能力（

我们要求

\[\int_L^R\dfrac{cx+d}{ax+b}\,\mathrm dx \]

记 \(u=ax+b\)，换元可得

\[\begin{aligned} &\int_L^R\dfrac{cx+d}{ax+b}\\ =&\int_L^R\dfrac{1}{a}·\dfrac{cx+d}{ax+b}·(ax+b)'\,\mathrm dx\\ =&\dfrac{1}{a}\int_L^R\dfrac{cx+d}{u}\,\mathrm du \end{aligned} \]

分母上以经出现了换元变量 \(u\)，考虑将分子也写成 \(u\) 的表达式：

\(cx+d=c·\dfrac{u-b}{a}+d=\dfrac{cu}{a}-(\dfrac{bc}{a}-d)\)

于是

\[\begin{aligned} &\dfrac{1}{a}\int_L^R\dfrac{cx+d}{u}\,\mathrm du\\ =&\dfrac{1}{a}\int_L^R\dfrac{\dfrac{cu}{a}-(\dfrac{bc}{a}-d)}{u}\,\mathrm du\\ =&\dfrac{1}{a}\int_L^R\dfrac{c}{a}-(\dfrac{bc}{a}-d)·\dfrac{1}{u}\,\mathrm du\\ =&\dfrac{1}{a}(\dfrac{c}{a}\int_L^R\,\mathrm du-(\dfrac{bc}{a}-d)\int_L^R\dfrac{1}{u}\,\mathrm du)\\ \end{aligned} \]

根据 \(\int_{a}^b\dfrac{1}{x}\,\mathrm dx=\ln|b|-\ln|a|\)

可得

\[\begin{aligned} &\dfrac{1}{a}(\dfrac{c}{a}\int_L^R\,\mathrm du-(\dfrac{bc}{a}-d)\int_L^R\dfrac{1}{u}\,\mathrm du)\\ =&g(R)-g(L) \end{aligned} \]

其中

\[g(x)=\dfrac{1}{a}(\dfrac{c}{a}(ax+b)-(\dfrac{bc}{a}-d)\ln|ax+b|) \]

当然要注意特判 \(a=0\) 的情况

\[\begin{aligned} &\int_L^R\dfrac{cx+d}{b}\,\mathrm dx\\ =&\dfrac{c}{b}\int_L^Rx\,\mathrm dx+\dfrac{d}{b}\int_L^R1\,\mathrm dx\\ =&\dfrac{c}{2b}(R^2-L^2)+\dfrac{d}{b}(R-L) \end{aligned} \]

P4526【模板】自适应辛普森法2

出 现 了 正 无 穷 ？ ！ ！ 反 常 积 分 ？ ！ ！

不过注意到一件事情，那就是对于 \(a<0\) 的情况，当 \(x\) 无限逼近 \(0\) 时函数值增长是很快的，于是我们猜测 \(a<0\) 是积分发散，有位大佬也给出了详细证明，然鹅我没看懂/ll。

当 \(a\ge 0\) 时候不难发现当 \(x\) 很大时 \(y\) 是非常非常接近 \(0\) 的，就按 \(a=50\) 来算，当 \(x=20\) 时 \(y=20^{-17.5}\approx 1.7\times 10^{-23}\)，这远远小于精度误差，因此我们可以把积分的上界直接当作 \(20\)，由于下界不能为 \(0\)，因此我们积分下界可以设为 \(\epsilon\)，按照上一题的套路来就行了。