辛普森积分
自适应辛普森积分是一种解决定积分求解问题的算法。
给出一个函数(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;
}