算法积分——自适应辛普森

手算是积分的最好办法。但是电脑不能手算！
对于离散的函数，可以直接累加输出量 * dt来完成。

C++

double f(double x){
	return //...
}
double dt = 1e-5;
double calc(double l,double r){
	double ans = 0.0;
	for(auto i=l;i<=r;i+=dt){
		ans += f(i)*dt;
	}
	return ans;
}

但是对于连续的函数，这种方法显然是不可行的。
于是考虑其他办法：

所以我们现在只需要考虑如何把函数分成多个二次函数即可。
那么怎么拟合呢？

\begin{aligned}&对于每个函数f(x)，它都能表示为由很多个\textbf{二次函数}组成的分段函数，又叫做\textbf{拟合}\\&对每一个二次函数，我们可以计算它的积分:\\&\int_{a}^{b} ax^{2}+ b x +c\, dx\\&=(\frac{a}{3}x^{3}+\frac{b}{2}x^2+cx)|^{b}_{a}\\&=\frac{b-a}{6}[f(a)+f(b)+4\cdot f(\frac{a+b}{2})]\end{aligned}

对整段，左半部分和右半部分用二次函数积分公式计算，结果记作 $A,L,R$ ，如果 $L+R \approx A$ 则认为这一段是与二次函数相似。

C++

#define mid ((l+r)/2)
double f(double x){
	return //...
}
double simpson(double l,double r){
	return (r-l)*(f(l)+f(r)+4*mid)/6;
}
double asr(double l,double r,double eps,double A){
	double L = simpson(l,mid),R = simpson(mid,r);
	if(fabs(L+R-A)<=15 * eps) return L + R + (L + R - A )/ 15;
	return asr(l,mid,eps/2,L)+asr(mid,r,eps/2,R)
}
int main(){
	return asr(l,r,1e-15,simpson(l,r));
}

例题

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

当然除了拟合成二次函数，一次函数/三次函数也是可行的。当然那就是其他的算法了。

算法积分——自适应辛普森

留言 / 评论