自适应辛普森法

@Pelom  October 11, 2021

自适应辛普森法

一、定积分

定积分是积分的一种,是函数$f(x)$在区间$[a,b]$上积分和的极限

定积分的几何意义是函数的曲线上$x$的一段区间与$x$轴围成的曲边梯形的带符号面积

表示为$$\int_a^b f(x)\mathrm{d}x$$
由定积分的定义,若函数$f(x)$在区间$[a,b]$上可积分,则有$n$等分的特殊分法$$\lim_{n \rightarrow +\infty} \sum_{i=1}^n f[a+\frac{i}{n}(b-a)]\frac{b-a}{n}=\int_a^b f(x)\mathrm{d}x$$

二、牛顿-莱布尼茨公式

牛顿-莱布尼兹公式(Newton-Leibniz formula),通常也被称为微积分基本定理,揭示了定积分与被积函数的原函数或者不定积分之间的联系

如果函数$f(x)$在区间$[a,b]$上连续,并且存在原函数$F(x)$
则$$\int_a^b f(x)\mathrm{d}x=F(b)-F(a)=F(x)|_a^b$$
原函数$F(x)$为满足$$F'(x)=f(x)$$

三、辛普森公式

辛普森(Simpson)公式是牛顿-科特斯公式当$n=2$时的情形,也称为三点公式。利用区间二等分的三个点来进行积分插值。其科特斯系数分别为$\frac{1}{6},\frac{4}{6},\frac{1}{6}$。

即用二次函数来拟合被积函数
设拟合后的二次函数为$$g(x)=Ax^2+Bx+C$$
则$g(x)$的原函数为$$G(x)=\frac{A}{3}x^3+\frac{B}{2}x^2+Cx$$

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

所以$$\int_a^b f(x)\mathrm{d}x \approx \frac{(b-a)(f(a)+4f(\frac{a+b}{2})+f(b))}{6}$$

代码:

inline double simpson(double l,double r){
    return (r-l)*(f(l)+4*f((l+r)/2)+f(r))/6; //f(x)为被积函数
}

四、自适应辛普森法

由于是用二次函数拟合原函数,自然会存在精度问题;而提高精度的方法是分为更多的小区间求和
根据给定的精度自动控制分区间的数量,这便是自适应辛普森法
算法由递归实现:若满足精度要求,则终止递归并返回;否则,二分区间递归求和

代码:

double asr(double l,double r,double res){
    double mid=(l+r)/2;
    double L=simpson(l,mid),R=simpson(mid,r);
    if(fabs(L+R-res)<=15*eps)
        return L+R+(L+R-res)/15;
    return asr(l,mid,L)+asr(mid,r,R);
}
inline double asr(double l,double r){
    return asr(l,r,simpson(l,r));
}

其中控制精度时出现的常数$15$由下图论证算出

注: 算法中的eps越小越精确,也更费时,尽量比题目要求精度低$4$个数量级


添加新评论