数值积分 定积分的定义 简单来说，函数 f(x) [l,r] \int_{l}^{r}f(x)\mathrm{d}x f(x) [l,r] x x x
很多情况下，我们需要高效，准确地求出一个积分的近似值。下面介绍的 辛普森法 ，就是这样一种求数值积分的方法。
辛普森法 这个方法的思想是将被积区间分为若干小段，每段套用二次函数的积分公式进行计算。
二次函数积分公式（辛普森公式） 对于一个二次函数 f(x)=Ax^2+Bx+C
\int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6} 推导过程： 对于一个二次函数 f(x)=Ax^2+Bx+C F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D
\begin{aligned} \int_l^r f(x) {\mathrm d}x &= F(r)-F(l) \\ &= \frac{a}{3}(r^3-l^3)+\frac{b}{2}(r^2-l^2)+c(r-l) \\ &=(r-l)(\frac{a}{3}(l^2+r^2+lr)+\frac{b}{2}(l+r)+c) \\ &=\frac{r-l}{6}(2al^2+2ar^2+2alr+3bl+3br+6c)\\ &=\frac{r-l}{6}((al^2+bl+c)+(ar^2+br+c)+4(a(\frac{l+r}{2})^2+b(\frac{l+r}{2})+c)) \\ &=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2})) \end{aligned} 根据这个辛普森公式，我们先介绍一种普通的辛普森积分法。
普通辛普森法 1743 年，这种方法发表于托马斯·辛普森的一篇论文中。
描述 给定一个自然数 n [l, r] 2n x
x_i = l + i h, ~~ i = 0 \ldots 2n, h = \frac {r-l} {2n}.
我们就可以计算每个小区间 [x_ {2i-2}, x_ {2i}] i = 1\ldots n
对于 [x_ {2i-2}, x_ {2i}] i = 1\ldots n (x_ {2i-2}, x_ {2i-1}, x_ {2i}) P(x) P(x)
\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}
将其分段求和即可得到如下结论：
\int_l ^ r f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}
误差 我们直接给出结论，普通辛普森法的误差为：
-\tfrac{1}{90} \left(\tfrac{r-l}{2}\right)^5 f^{(4)}(\xi) 其中 \xi [l,r]
实现 1
2
3
4
5
6
7
8
9
10
11
12
13 // C++ Version
const int N = 1000 * 1000 ; double simpson_integration ( double a , double b ) { double h = ( b - a ) / N ; double s = f ( a ) + f ( b ); for ( int i = 1 ; i <= N - 1 ; ++ i ) { double x = a + h * i ; s += f ( x ) * (( i & 1 ) ? 4 : 2 ); } s *= h / 3 ; return s ; }
1
2
3
4
5
6
7
8
9
10
11
12
13
14 # Python Version
N = 1000 * 1000
def simpson_integration ( a , b ):
h = ( b - a ) / N
s = f ( a ) + f ( b )
for i in range ( 1 , N ):
x = a + h * i
if i & 1 :
s = s + f ( x ) * 4
else :
s = s + f ( x ) * 2
s = s * ( h / 3 )
return s
自适应辛普森法 普通的方法为保证精度在时间方面无疑会受到 n
现在唯一的问题就是如何进行分段。如果段数少了计算误差就大，段数多了时间效率又会低。我们需要找到一个准确度和效率的平衡点。
我们这样考虑：假如有一段图像已经很接近二次函数的话，直接带入公式求积分，得到的值精度就很高了，不需要再继续分割这一段了。
于是我们有了这样一种分割方法：每次判断当前段和二次函数的相似程度，如果足够相似的话就直接代入公式计算，否则将当前段分割成左右两段递归求解。
现在就剩下一个问题了：如果判断每一段和二次函数是否相似？
我们把当前段直接代入公式求积分，再将当前段从中点分割成两段，把这两段再直接代入公式求积分。如果当前段的积分和分割成两段后的积分之和相差很小的话，就可以认为当前段和二次函数很相似了，不用再递归分割了。
上面就是自适应辛普森法的思想。在分治判断的时候，除了判断精度是否正确，一般还要强制执行最少的迭代次数。
参考代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 // C++ Version
double simpson ( double l , double r ) { double mid = ( l + r ) / 2 ; return ( r - l ) * ( f ( l ) + 4 * f ( mid ) + f ( r )) / 6 ; // 辛普森公式
} double asr ( double l , double r , double eqs , double ans , int step ) { double mid = ( l + r ) / 2 ; double fl = simpson ( l , mid ), fr = simpson ( mid , r ); if ( abs ( fl + fr - ans ) <= 15 * eqs && step < 0 ) return fl + fr + ( fl + fr - ans ) / 15 ; // 足够相似的话就直接返回
return asr ( l , mid , eqs / 2 , fl , step - 1 ) + asr ( mid , r , eqs / 2 , fr , step - 1 ); // 否则分割成两段递归求解
} double calc ( double l , double r , double eps ) { return asr ( l , r , eps , simpson ( l , r ), 12 ); }
1
2
3
4
5
6
7
8
9
10
11
12
13 # Python Version
def simpson ( l , r ):
mid = ( l + r ) / 2
return ( r - l ) * ( f ( l ) + 4 * f ( mid ) + f ( r )) / 6 # 辛普森公式
def asr ( l , r , eqs , ans , step ):
mid = ( l + r ) / 2
fl = simpson ( l , mid ); fr = simpson ( mid , r )
if abs ( fl + fr - ans ) <= 15 * eqs and step < 0 :
return fl + fr + ( fl + fr - ans ) / 15 # 足够相似的话就直接返回
return asr ( l , mid , eqs / 2 , fl , step - 1 ) + \
asr ( mid , r , eqs / 2 , fr , step - 1 ) # 否则分割成两段递归求解
def calc ( l , r , eps ):
return asr ( l , r , eps , simpson ( l , r ), 12 )
习题 build 本页面最近更新：更新历史 edit 发现错误？想一起完善？ 在 GitHub 上编辑此页！ people 本页面贡献者：OI-wiki copyright 本页面的全部内容在 CC BY-SA 4.0 和 SATA