利用辛普森公式求积分（Integration by Simpson’s formula）
我们将计算定积分的值 ∫ a b f ( x ) d x \int_a ^ b f (x) dx ∫abf(x)dx 。
这里描述的解决方案发表在 1743 年Thomas Simpson的一篇论文中。
辛普森公式
令 n n n 是某自然数。我们将积分段 [ a , b ] [a, b] [a,b] 分成 2 n 2n 2n 个相等的部分：
x i = a + i h , i = 0 … 2 n , h = b − a 2 n . x_i = a + i h, ~~ i = 0 \ldots 2n,~~h = \frac {b-a} {2n}. xi=a+ih, i=0…2n, h=2nb−a.
现在我们分别计算每个段上的积分 [ x 2 i − 2 , x 2 i ] [x_ {2i-2}, x_ {2i}] [x2i−2,x2i] ， i = 1 … n i = 1 \ldots n i=1…n ，然后添加所有值。
所以，假设我们考虑下一个线段 [ x 2 i − 2 , x 2 i ] , i = 1 … n [x_ {2i-2}, x_ {2i}], i = 1 \ldots n [x2i−2,x2i],i=1…n 。用经过3个点 ( x 2 i − 2 , x 2 i − 1 , x 2 i ) (x_ {2i-2}, x_ {2i-1}, x_ {2i}) (x2i−2,x2i−1,x2i) 的抛物线 P ( x ) P(x) P(x) 代替其上的函数 f ( x ) f(x) f(x) 。这样的抛物线一直存在，并且是唯一的；它可以通过分析找到。例如，我们可以使用拉格朗日多项式插值来构造它。剩下要做的唯一事情就是整合这个多项式。如果你对一个一般的函数 f f f 这样做，你会得到一个非常简单的表达式：
∫ x 2 i − 2 x 2 i f ( x ) d x ≈ ∫ x 2 i − 2 x 2 i P ( x ) d x = ( f ( x 2 i − 2 ) + 4 f ( x 2 i − 1 ) + ( f ( x 2 i ) ) h 3 \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} ∫x2i−2x2if(x) dx≈∫x2i−2x2iP(x) dx=(f(x2i−2)+4f(x2i−1)+(f(x2i))3h
在所有段上添加这些值，我们得到最终的辛普森公式：
∫ a b f ( x ) d x ≈ ( f ( x 0 ) + 4 f ( x 1 ) + 2 f ( x 2 ) + 4 f ( x 3 ) + 2 f ( x 4 ) + … + 4 f ( x 2 N − 1 ) + f ( x 2 N ) ) h 3 \int_a ^ b 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} ∫abf(x)dx≈(f(x0)+4f(x1)+2f(x2)+4f(x3)+2f(x4)+…+4f(x2N−1)+f(x2N))3h
误差
用辛普森公式逼近积分的误差是
− 1 90 ( b − a 2 ) 5 f ( 4 ) ( ξ ) -\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi) −901(2b−a)5f(4)(ξ)
其中 ξ \xi ξ 是介于 a a a 和 b b b 之间的数。
实现
代码中， f ( x ) f(x) f(x) 代表用户定义的函数。
const int N = 1000 * 1000; // number of steps (already multiplied by 2) double simpson_integration(double a, double b){
double h = (b - a) / N; double s = f(a) + f(b); // a = x_0 and b = x_2n for (int i = 1; i <= N - 1; ++i) {
// Refer to final Simpson's formula double x = a + h * i; s += f(x) * ((i & 1) ? 4 : 2); } s *= h / 3; return s; } 发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/216582.html原文链接：https://javaforall.net
