自适应辛普森法初探

求函数定积分用，以两道咕咕上的题目为例。

题一【p4525】

直接套用“牛顿莱布尼茨”公式，被积函数
\[ f(x) =\frac{cx+d}{ax+b} =\frac{c}{a}(1+\frac{\frac{ad}{c}-b}{ax+b})\ \]
情况一，\(ac\not=0​\)
\[ \int f(x)dx =\frac{c}{a}(\int dx+(\frac{ad}{c}-b)\int\frac{1}{ax+b}dx)+c\\ =\frac{c}{a}(x+(\frac{ad}{c}-b)(\frac{1}{a}\ln|ax+b|))+c\\ =\frac{c}{a}x+(d-\frac{bc}{a})\frac{1}{a}\ln|ax+b|+c\\ \]
情况二，\(a=0​\)
\[ \int f(x)dx =\frac{1}{b}\int(cx+d)dx=\frac{1}{b}(\frac{cx^2}{2}+dx)+c \]
情况三，\(c=0​\)
\[ \int f(x)dx=d\int\frac{1}{ax+b}dx=\frac{d}{a}\ln|ax+b|+c \]
情况二和情况三都能计算\(a=c=0\)的答案。

#include <bits/stdc++.h>
using namespace std;
double a,b,c,d,l,r;

double f1(double x) {
    return c/a*x+(d-b*c/a)/a*log(fabs(a*x+b));
}
double f2(double x) {
    return (c*x*x/2+d*x)/b;
}
double f3(double x) {
    return d*log(fabs(a*x+b))/a;
}


int main() {
    cin>>a>>b>>c>>d>>l>>r;
    if(a*c!=0) printf("%.6f\n",f1(r)-f1(l));
    else if(!a) printf("%.6f\n",f2(r)-f2(l));
    else printf("%.6f\n",f3(r)-f3(l));
    return 0;
}

题二【p4526】

simpson公式的推导：用\(g(x)=ax^2+bx+c\)来拟合\(f(x)\)，则有
\[ \int f(x)dx\approx\int g(x)dx=\int(ax^2+bx+c)dx=\frac{a}{3}x^3+\frac{b}{2}x^2+cx+c \]
由“牛顿莱布尼茨”公式
\[ \int_{a}^b f(x)dx\approx\frac{a}{3}(b^3-a^3)+\frac{b}{2}(b^2-a^2)+c(b-a)\\ =\frac{a}{3}(b-a)(b^2+ab+a^2)+\frac{b}{2}(b-a)(b+a)+c(b-a)\\ =\frac{b-a}{6}(2ab^2+2aab+2aa^2+3bb+3ba+6c)\\ =\frac{b-a}{6}(aa^2+ba+c+ab^2+bb+c+ab^2+a2ab+aa^2+bb+ba+4c)\\ =\frac{b-a}{6}(f(a)+f(b)+4f(\frac{a+b}{2})) \]
然后是自适应辛普森算法（asr），每次用simpson公式算出当前区间的左、右部分之和，如果和整个区间的公式值相差无几，则认为此范围内精度足够，退出迭代。

当\(a<0\)时，\(\lim_{x\to0}f(x)=-\infty\)，题目函数发散。

#include <bits/stdc++.h>
using namespace std;
const double eps=1e-9;

double a;
double f(double x) {
    return pow(x,a/x-x);
}
double simpson(double a,double b) {
    return (b-a)*(f(b)+f(a)+4*f((a+b)/2))/6;
}
double ars(double l,double r,double tot) {
    double mid=(l+r)/2;
    double ls=simpson(l,mid);
    double rs=simpson(mid,r);
    if(fabs(ls+rs-tot)<eps) return ls+rs;
    return ars(l,mid,ls)+ars(mid,r,rs);
}
double ars(double l,double r) {
    return ars(l,r,simpson(l,r));
}

int main() {
    scanf("%lf",&a);
    if(a<0) puts("orz");
    else printf("%.5f\n",ars(eps,30));
    return 0;
}