自适应辛普森积分

辛普森积分

自适应辛普森积分是一种解决定积分求解问题的算法。

给出一个函数(f(x))，求：

[int _l^rf(x){ m{d}}x ]

我们考虑用一条抛物线来近似这个函数，设(g(x)=ax^2+bx+c)。

那么可得：

[egin{align} int_l^rf(x){ m{d}}x&approx int_l^rg(x){ m{d}}x\ &=frac{1}{3}a(r^3-l^3)+frac{1}{2}b(r^2-l^2)+c(r-l)\ &=frac{(r-l)(al^2+bl+c+ar^2+br+c+a(l+r)^2+2b(l+r)+4c}{6}\ &=frac{(r-l)(f(l)+f(r)+4f(frac{l+r}{2}))}{6} end{align} ]

那么这个玩意就叫(simpson)公式。

自适应辛普森积分

上面这个玩意显然精度不够，甚至可以说根本就不对。

那么怎么解决这个问题呢，我们可以考虑模拟微分的过程，把定义域切成一小段一小段，然后在近似，这样精度就有了。

但是显然这样是很慢的，所以我们有了一种自适应的想法，即若当前区间精度够了就退出，否则二分然后递归计算左右两个区间，具体来说，代码应该这么写：

double _int (double l,double r) {return (r-l)*(f(l)+f(r)+4.0*f((l+r)/2.0))/6.0;}
double calc(double l,double r,double ans) {
	double mid=(l+r)/2.0,L=_int(l,mid),R=_int(mid,r);
	if(fabs(L+R-res)<eps) return res;
	else return calc(l,mid,L)+calc(mid,r,R);
}

但是这样精度可能还是不够，然后有一种更玄学的近似，具体证明我也不会...

代码具体长这样：

double calc(double l,double r,double res,double Eps) {
	double mid=(l+r)/2.0,L=_int(l,mid),R=_int(mid,r);
	if(fabs(L+R-res)<Eps*15.0) return L+R+(L+R-res)/15.0;
	else return calc(l,mid,L,Eps*0.5)+calc(mid,r,R,Eps*0.5);
}

例题

题目链接：luogu P4525。

这题就直接套公式就好了。

#include<bits/stdc++.h>
using namespace std;
 
void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double 

const int maxn = 2e5+10;
const lf eps = 1e-9;

lf a,b,c,d;

lf f(lf x) {return (c*x+d)/(a*x+b);}

lf integral(lf L,lf R) {return (R-L)*(f(L)+f(R)+4.0*f((L+R)/2.0))/6.0;}

lf calc(lf L,lf R,lf ans,lf Eps) {
	lf mid=(L+R)/2.0;
	lf l=integral(L,mid),r=integral(mid,R);
	if(fabs(l+r-ans)<Eps*15.0) return l+r+(l+r-ans)/15.0;
	else return calc(L,mid,l,Eps*0.5)+calc(mid,R,r,Eps*0.5);
}

int main() {
	lf L,R;
	scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
	printf("%.6lf
",calc(L,R,integral(L,R),eps));
	return 0;
}

不过这个积分很容易积出来：

[egin{align} intfrac{ax+b}{cx+d}{ m{d}}x&=intfrac{c}{a}-frac{d-frac{bc}{a}}{ax+b}{ m{d}}x\ &=frac{cx}{a}-frac{1}{a}(d-frac{bc}{a})ln|ax+b|+C end{align} ]

注意到中间有一个很简单的换元：

[int frac{1}{ax+b}{ m{d}}x=frac{1}{a}intfrac{1}{ax+b}{ m{d}}(ax+b)=frac{1}{a}ln|ax+b| ]

---

题目链接：luogu P4526。

注意到(a<0)时，(lim_{x o 0}f(x)=+infty)，此时积分发散。

否则积分收敛，我也不是特别清楚为什么

那么直接套公式就好了。

#include<bits/stdc++.h>
using namespace std;
 
void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double 

const int maxn = 2e5+10;
const lf eps = 1e-8;

lf a;

lf f(lf x) {return pow(x,a/x-x);}
lf _int(lf l,lf r) {return (r-l)*(f(l)+f(r)+4.0*f((l+r)/2.0))/6.0;}
lf calc(lf l,lf r,lf res,lf Eps) {
	lf mid=(l+r)/2.0,L=_int(l,mid),R=_int(mid,r);
	if(fabs(L+R-res)<eps*15.0) return L+R+(L+R-res)/15.0;
	else return calc(l,mid,L,Eps*0.5)+calc(mid,r,R,Eps*0.5);
}

int main() {
	scanf("%lf",&a);
	if(a<0) puts("orz");
	else printf("%.5lf
",calc(eps,20.0,_int(eps,20.0),eps));
	return 0;
}