辛普森积分
参考资料
知乎_数值积分漫谈 （推荐阅读）
前置
牛顿-莱布尼茨公式（积分基本公式）
普通的Simpson积分
一般的，我们用小的矩阵或梯形来近似的代表一小段的积分值，如
当然， \(a\) 和 \(b\) 之间的距离应当足够小。
这实际上是用过点 \((a,f(a))\) 和点 \((b,f(b))\) 的直线代替了这一段的函数。
为了增加效率并减小误差，Simpson积分用二次函数代替一段函数。
设替代的函数为 \(g(x)=Ax^2+Bx+C\) ，我们需要知道这三个值： \(g(a)\) ， \(g(b)\) 和 \(g(\frac{a+b}{2})\) 。
那么根据牛顿-莱布尼茨公式，
所以我们只要在将一个区间分成足够多块，就能求出精度足够高的值。
自适应Simpson积分（Adaptive Simpson's method）
要达到较高的精度往往会耗费很多的时间。而如果函数的某一部分比较平滑，那可以将这个区间少分点段；同样的，将不太平滑的区间多分几段。
我们可以二分计算区间积分，同时通过期望容差来控制二分的终止。
一般的，如果满足 $|S(a,m)+S(m,b)-S(a,b)|<15\epsilon $ ，就终止二分并返回 \(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+k\frac{b-a}{n}\) ，那么
其中
当 \(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\n", 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\n" );
else printf( "%.5lf\n", Integral( 4e-7, Inf, 4e-7, Simpson( 4e-7, Inf ) ) );
return;
}
