自适应辛普森法

自适应辛普森法通过 $Simpson$ 公式，用二次函数来拟合，实现时用二分递归来自动控制区间分割的大小，既保证精度，又保证速度

$Simpson$ 公式推导

\int_{a}^{b} f (x) d x

$\int_a^bf(x)dx$

\approx \int_{a}^{b} A x^{2} + B x + C

$\approx\int_a^bAx^2+Bx+C$

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

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

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

$=\frac{(b-a)}{6}(2Ab^2+2Aab+2Aa^2+3Bb+3Ba+6C)$

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

$=\frac{(b-a)}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C)$

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

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

求 $\displaystyle{\int_l^r\frac{cx+d}{ax+b}\mathrm{d}x}$

$code:$

double f(double x)
{
    return (c*x+d)/(a*x+b);
}
double simpson(double l,double r)
{
    double mid=(l+r)/2;
    return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
double solve(double l,double r,double eps)
{
    double mid=(l+r)/2;
    double a=simpson(l,mid),b=simpson(mid,r),m=simpson(l,r);
    if(fabs(a+b-m)<=15*eps) return a+b+(a+b-m)/15;
    return solve(l,mid,eps/2)+solve(mid,r,eps/2);
}

月下柠檬树：求一个图形的面积

$code:$

double f(double x)
{
    double y=0;
    for(int i=1;i<n;++i)
        if(x>=h[i]-r[i]&&x<=h[i]+r[i])
            y=max(y,sqrt(r[i]*r[i]-(h[i]-x)*(h[i]-x)));
    for(int i=1;i<n;++i)
        if(x>=p[i].l&&x<=p[i].r)
            y=max(y,p[i].k*x+p[i].b);
    return y;
}
double simpson(double l,double r)
{
    double mid=(l+r)/2;
    return (r-l)*(f(l)+f(r)+4*f(mid))/6;
}
double solve(double l,double r,double eps)
{
    double mid=(l+r)/2;
    double a=simpson(l,mid),b=simpson(mid,r),m=simpson(l,r);
    if(fabs(a+b-m)<=15*eps) return a+b+(a+b-m)/15;
    return solve(l,mid,eps/2)+solve(mid,r,eps/2);
}