数值分析——自适应辛普森积分
自适应辛普森积分法是利用二次函数对被积函数进行插值的一种模拟积分法,常用于估计积分的数值,利用递归自动适应所需精度。
公式推导
基本思想公式:
∫ a b f ( x ) d x ≈ ∫ a b ( A x 2 + B x + C ) d x \int_{a}^{b} f(x) dx \approx \int_{a}^{b} (Ax^2+Bx+C)dx ∫abf(x)dx≈∫ab(Ax2+Bx+C)dx
然后我们对这个二次函数积分进行求解:
∫ a b ( A x 2 + B x + C ) d x = ( b − a ) 6 ( A a 2 + B a + C + A b 2 + B b + C + 4 ( A ( a + b 2 ) 2 + B ( a + b 2 ) + C ) ) ) \int_{a}^{b} (Ax^2+Bx+C)dx = \frac{(b-a)}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4(A(\frac{a+b}2)^2+B(\frac{a+b}2)+C))) ∫ab(Ax2+Bx+C)dx=6(b−a)(Aa2+Ba+C+Ab2+Bb+C+4(A(2a+b)2+B(2a+b)+C)))
其中,我们用:
f ( a ) = A a 2 + B a + C f ( b ) = A b 2 + B b + C 4 f ( a + b 2 ) = 4 ( A ( a + b 2 ) 2 + B ( a + b 2 ) + C ) ) \begin{matrix} f(a) & = & Aa^2+Ba+C \\ f(b) & = & Ab^2+Bb+C \\ 4f(\frac{a + b}{2}) & = & 4(A(\frac{a+b}2)^2+B(\frac{a+b}2)+C)) \end{matrix} f(a)f(b)4f(2a+b)===Aa2+Ba+CAb2+Bb+C4(A(2a+b)2+B(2a+b)+C))
去替换原式得到最终的二次自适应辛普森积分公式:
∫ a b f ( x ) d x ≈ ∫ a b ( A x 2 + B x + C ) d x = ( b − a ) 6 ( f ( a ) + f ( b ) + 4 f ( a + b 2 ) ) \int_{a}^{b} f(x) dx \approx \int_{a}^{b} (Ax^2+Bx+C)dx = \frac{(b-a)}{6} (f(a) + f(b) + 4f(\frac{a + b}{2})) ∫abf(x)dx≈∫ab(Ax2+Bx+C)dx=6(b−a)(f(a)+f(b)+4f(2a+b))
使用
首先我们定义我们的被积函数:
double f(double x)
{
return (c * x + d) / (a * x + b);
}
其次编写我的辛普森积分公式:
double simpson(double l, double r)
{
return (r - l) / 6 * (f(r) + f(l) + 4 * f((l + r) / 2));
}
然后是自适应精度进行递归计算:
double asr(double L, double R, double ans)
{
// 取中点
double mid = (L + R) / 2;
// 分别计算中点两边的积分数值
double l = simpson(L, mid);
double r = simpson(mid, R);
// 如果和参考答案的精度适应,那么直接返回积分值
if (abs(l + r - ans) <= esp)
return l + r;
// 否则将l,r作为参考答案,递归计算两边的积分数值
return asr(L, mid, l) + asr(mid, R, r);
}
说明
- 如果积分收敛,初始参考答案可理论可设置为任意值不影响结果,可以用中值作为积分值或者直接带入simpson作为积分值。
- 该方法的时间复杂度为 O ( n log ( n / e s p ) ) O(n \log (n/esp)) O(nlog(n/esp)), n n n为积分区间的长度。