辛普森积分

参考资料
前置
普通的Simpson积分
自适应Simpson积分（Adaptive Simpson's method）
牛顿-科斯特公式
板子

[ ewcommand{d}{mathrm{d}\,} ]

参考资料

百度百科_牛顿-莱布尼茨公式

知乎_数值积分漫谈（推荐阅读）

前置

牛顿-莱布尼茨公式（积分基本公式）

[int_a^bf(x) d x=F(b)-F(a)=F(x)|_a^b ]

普通的Simpson积分

一般的，我们用小的矩阵或梯形来近似的代表一小段的积分值，如

[int_a^bf(x) d x approx frac{(f(a)+f(b))(b-a)}{2} ]

当然， (a) 和 (b) 之间的距离应当足够小。

这实际上是用过点 ((a,f(a))) 和点 ((b,f(b))) 的直线代替了这一段的函数。

为了增加效率并减小误差，Simpson积分用二次函数代替一段函数。

设替代的函数为 (g(x)=Ax^2+Bx+C) ，我们需要知道这三个值： (g(a)) ， (g(b)) 和 (g(frac{a+b}{2})) 。

那么根据牛顿-莱布尼茨公式，

[egin{aligned} int_a^b f(x) d x &approx int_a^b g(x) d x\ &=G(b)-G(a)\ &=frac{A}{3}(a^3-b^3)+frac{A}{2}(a^2-b^2)+C(a-b)\ &=frac{A}{3}(a-b)(a^2+ab+b^2)+frac{B}{2}(a-b)(a+b)+C(a-b)\ &=frac{a-b}{6}(2Aa^2+2Aab+2Ab^2+3Ba+3Bb+6C)\ &=frac{a-b}{6}[(Aa^2+Ba+C)+(Ab^2+Bb+C)+4(Afrac{a^2+b^2+2ab}{4}+Bfrac{a+b}{2}+C)]\ &=frac{a-b}{6}(g(a)+g(b)+g(a+b))\ &=frac{a-b}{6}(g(a)+g(b)+4g(frac{a+b}{2})) end{aligned} ]

所以我们只要在将一个区间分成足够多块，就能求出精度足够高的值。

自适应Simpson积分（Adaptive Simpson's method）

要达到较高的精度往往会耗费很多的时间。而如果函数的某一部分比较平滑，那可以将这个区间少分点段；同样的，将不太平滑的区间多分几段。

我们可以二分计算区间积分，同时通过期望容差来控制二分的终止。

一般的，如果满足 $|S(a,m)+S(m,b)-S(a,b)|<15epsilon $ ，就终止二分并返回 (S(a,m)+S(m,b)+frac{S(a,m)+S(m,b)-S(a,b)}{15})。

其中 (m=frac{a+b}{2}) ， (epsilon) 为期望容差。

所以可以这样实现：

double F(double x);

double Simpson( double L, double R ) {
    double Mid = (L + R) / 2.0;
    return ( F( L ) + 4 * F( Mid ) + F( R ) ) * ( R - L ) / 6;
}

double Integral( double L, double R, double Eps ) {
    double Mid = (L + R) / 2.0;
    double ST = Simpson(L, R), SL = Simpson(L, Mid), SR = Simpson(Mid, R);
    if( std::fabs( SL + SR - ST ) <= 15 * Eps )  return SL + SR + (SL + SR - ST) / 15;
    return Integral( L, Mid, Eps / 2 ) + Integral( Mid, R, Eps / 2 );
}

一般，OI & ACM中要用到的到这里就差不多了。下面的内容仅做介绍~~毕竟我也不会~~。

牛顿-科斯特公式

假设 (I=[a,b]) ， (x_k=a+kfrac{b-a}{n}) ，那么

[int_If(t)d tapprox I_{appr}(x_1,x_2,cdots,x_n)=(b-a)sumlimits_{k=0}^nC_k^{(n)}f(x_k) ]

其中

[C_k^{(n)}=frac{1}{n}int_0^nprodlimits_{k eq j}frac{t-j}{k-j}d t=frac{(-1)^{n-k}}{ncdot k!(n-k)!}int_0^nprod_{k eq j}(t-j)d t ]

当 (n=1) 时，这个公式就是上面提到的梯形近似方法；而当 (n=2) 时，就是辛普森积分。

当 (n=3) 时，这个公式为 Simpson's 3/8 rule ，用更大的常数（在平滑函数下）换更高的精度~~手动狗头~~。

板子

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

#include <cstdio>
#include <algorithm>
#include <cmath>

inline void Init();
inline void Calc();
int main() { return Init(), Calc(), 0; }

double a, b, c, d, L, R;
inline void Init() { scanf( "%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &L, &R ); return; }
inline double f( double x ) { return ( c * x + d ) / ( a * x + b ); }
inline double Simpson( double l, double r ) { return ( r - l ) / 6 * ( f( l ) + f( r ) + 4 * f( ( l + r ) / 2 ) ); }
double Integral( double l, double r, double Eps, double Total ) {
	double Mid = ( l + r ) / 2;
	double L = Simpson( l, Mid ), R = Simpson( Mid, r );
	if( std::fabs( L + R - Total ) < 15 * Eps ) return L + R + ( L + R - Total ) / 15;
	return Integral( l, Mid, Eps / 2, L ) + Integral( Mid, r, Eps / 2, R );
}
inline void Calc() { printf( "%.6lf
", Integral( L, R, 4e-8, Simpson( L, R ) ) ); return; }

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

(f(0)) 是没有用的，先枪毙掉。首先，当 (a<0) 时，若 (x) 趋向于 (0) ，那么 (f(x)) 趋向于无穷，积分发散。

考虑 (a>0) ，发现函数收敛的很快。所以如果上界取得太大，会导致答案为 (0) 。所以直观感受一下，上界应该取得小一点。就Ok了。

#include <cstdio>
#include <algorithm>
#include <cmath>

inline void Init();
inline void Calc();
int main() { return Init(), Calc(), 0; }

const double Inf = 20;
double a;
inline void Init() { scanf( "%lf", &a ); return; }
inline double f( double x ) { return std::pow( x, a / x - x ); }
inline double Simpson( double l, double r ) { return ( r - l ) / 6 * ( f( l ) + f( r ) + 4 * f( ( l + r ) / 2 ) ); }
double Integral( double l, double r, double Eps, double Total ) {
	double Mid = ( l + r ) / 2;
	double L = Simpson( l, Mid ), R = Simpson( Mid, r );
	if( std::fabs( L + R - Total ) < 15 * Eps ) return L + R + ( L + R - Total ) / 15;
	return Integral( l, Mid, Eps / 2, L ) + Integral( Mid, r, Eps / 2, R );
}
inline void Calc() { 
	if( a < 0 ) printf( "orz
" ); 
	else printf( "%.5lf
", Integral( 4e-7, Inf, 4e-7, Simpson( 4e-7, Inf ) ) );
	return; 
}