自适应辛普森积分
辛普森积分
自适应辛普森积分是一种解决定积分求解问题的算法。
给出一个函数\(f(x)\)，求：
\[\int _l^rf(x){\rm{d}}x
\]
我们考虑用一条抛物线来近似这个函数，设\(g(x)=ax^2+bx+c\)。
那么可得：
\[\begin{align}
\int_l^rf(x){\rm{d}}x&\approx \int_l^rg(x){\rm{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('\n');}
#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\n",calc(L,R,integral(L,R),eps));
return 0;
}
不过这个积分很容易积出来：
\[\begin{align}
\int\frac{ax+b}{cx+d}{\rm{d}}x&=\int\frac{c}{a}-\frac{d-\frac{bc}{a}}{ax+b}{\rm{d}}x\\
&=\frac{cx}{a}-\frac{1}{a}(d-\frac{bc}{a})\ln|ax+b|+C
\end{align}
\]
注意到中间有一个很简单的换元：
\[\int \frac{1}{ax+b}{\rm{d}}x=\frac{1}{a}\int\frac{1}{ax+b}{\rm{d}}(ax+b)=\frac{1}{a}\ln|ax+b|
\]
题目链接：luogu P4526。
注意到\(a<0\)时，\(\lim_{x\to 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('\n');}
#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\n",calc(eps,20.0,_int(eps,20.0),eps));
return 0;
}
