自适应辛普森法初探

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

题一【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;
}