自适应辛普森积分法

今天谈一谈自适应辛普森积分法

积分定义不多说了

那辛普森积分法是干什么的呢？

我们想求(displaystyle int _{a}^{b}displaystyle f(x)dx)

首先假设你现在并不会直接求积分,我们只能求个近似值

那求近似值方法就多了,我们可以近似成矩形和,梯形和

但是这些不是时间复杂度太高,就是精度太低

我们想为什么用矩形,梯形拟合曲线时存在较大误差,显然用直线拟合曲线不是很靠谱

我们就想到使用曲线拟合曲线

为了是表达简单我们选择的曲线是二次函数

当然这也不是最好的方法,但他在实际运行效果还是不错的

更好的拟合方式是牛顿-科特斯公式什么的，由于我也不会,暂时不说了

设(f(x)approx Ax^2+Bx+C)

(=displaystyle int _{a}^{b}f(x)dxapprox displaystyle int _{a}^{b}Ax^2+Bx+C)

(=displaystyle int _{a}^{b}Ax^2+displaystyle int _{a}^{b}Bx+displaystyle int _{a}^{b}C)

(=displaystyle frac{A}{3}x^3+frac{B}{2}x^2+Cxdisplaystyle|_{a}^{b})

(=frac{(a-b)}{6}(2{A}(a^2+ab+b^2)+3{B}(a+b)+6C))

(=frac{(a-b)}{6}(({A}a^2+{B}a+C)+({A}b^2+{B}b+C)+4((frac{a+b}{2})^2A+(frac{a+b}{2})B+C)))

(=frac{(a-b)}{6}(f(a)+f(b)+4*f(frac{a+b}{2})))

于是得出最终结论(displaystyle int _{a}^{b}f(x)dx) (approx) (frac{(a-b)}{6}(f(a)+f(b)+4*f(frac{a+b}{2})))

当然光这么算精度差太多

我们边从中点将([a,b])分成([a,frac{a+b}{2}]),([frac{a+b}{2},b])

类似分治

当精度够了,也就是说(displaystyle int _{a}^{frac{a+b}{2}}f(x)dx)+(displaystyle int ^{b}_{frac{a+b}{2}}f(x)dx-displaystyle int _{a}^{b}f(x)dx<eps)
时便停止递归
(Code)

//f(x) 为原函数
const db eps=1.0e-8;
#define hps eps*0.5
inline db gral(db l,db r)
{
    return (r-l)/6.0*(f(l)+f(r)+4*f((l+r)/2.0));
}
inline db sim(db l,db r,db v,db eps)
{
    db mid=(l+r)/2.0,t;
    db lv=gral(l,mid),rv=gral(mid,r);
    if(fabs(t=lv+rv-v)<eps*15) return lv+rv+t/15;
    return sim(l,mid,lv,hps)+sim(mid,r,rv,hps);
}

自适应辛普森积分也算说完了