c/c++语言开发共享自适应辛普森法初探

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

题一【p4525】

直接套用“牛顿莱布尼茨”公式，被积函数
[ f(x) =frac{cx+d}{ax+b} =frac{c}{a}(1+frac{frac{ad}{c}-b}{ax+b}) ]
情况一，(acnot=0​)
[ int f(x)dx =frac{c}{a}(int dx+(frac{ad}{c}-b)intfrac{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=dintfrac{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("%.6fn",f1(r)-f1(l));     else if(!a) printf("%.6fn",f2(r)-f2(l));     else printf("%.6fn",f3(r)-f3(l));     return 0; }

题二【p4526】

simpson公式的推导：用(g(x)=ax^2+bx+c)来拟合(f(x))，则有
[ int f(x)dxapproxint 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)dxapproxfrac{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_{xto0}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("%.5fn",ars(eps,30));     return 0; }