我发现我纯属闲着没事干。
这玩意拿来求定积分。
辛普森公式
对于二次函数 \(f(x)=ax^2+bx+c\)，有：
\[\int_l^rf(x)\text dx=\frac{(r-l)(f(l)+f(r)+4f(\frac{l+r}2))}6 \]证明大力拆。dirty work。
普通辛普森法
给个自然数 \(n\)，把区间 \([l,r]\) 分成 \(2n\) 个等长区间，每个区间用辛普森公式拟合积分值，然后加起来。
oi-wiki 说这玩意误差是
\[-\frac 1{90}(\frac{r-l}2)^5f^(4)(\xi) \]其中 \(xi\) 是 \([l,r]\) 中的某个值。当然我们不用他。
自适应辛普森法
这玩意可以平衡时间和精度的限制。我们发现问题在怎么分段。
有一个简单明了的方法：如果这一段的图象和二次函数足够相似就停。
怎么判断？每次当前段带入公式求积分，然后从中点砍两半分别求积分。如果差值在精度误差内就停。
当然一般我们还要有一个最小迭代次数，迭代次数到了才停。
代码
double simpson(double l,double r){
double mid=(l+r)/2;
return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans,int dep){
double mid=(l+r)/2;
double fl=simpson(l,mid),fr=simpson(mid,r);
if(abs(fl+fr-ans)<=15*eps&&dep<0)return fl+fr+(fl+fr-ans)/15;
return asr(l,mid,eps/2,fl,dep-1)+asr(mid,r,eps/2,fr,dep-1);
}
double calc(double l,double r,double eps){
return asr(l,r,eps,simpson(l,r),12);
}
模板 1：
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
double a,b,c,d,L,R;
double f(double x){
return (c*x+d)/(a*x+b);
}
double simpson(double l,double r){
double mid=(l+r)/2;
return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans){
double mid=(l+r)/2;
double fl=simpson(l,mid),fr=simpson(mid,r);
if(abs(fl+fr-ans)<=15*eps)return fl+fr+(fl+fr-ans)/15;
return asr(l,mid,eps/2,fl)+asr(mid,r,eps/2,fr);
}
double calc(double l,double r,double eps){
return asr(l,r,eps,simpson(l,r));
}
int main(){
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
double ans=calc(L,R,1e-6);
printf("%.6lf\n",ans);
return 0;
}
模板 2：
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
double a;
double f(double x){
return pow(x,a/x-x);
}
double simpson(double l,double r){
double mid=(l+r)/2;
return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans,int dep){
double mid=(l+r)/2;
double fl=simpson(l,mid),fr=simpson(mid,r);
if(abs(fl+fr-ans)<=15*eps&&dep<0)return fl+fr+(fl+fr-ans)/15;
return asr(l,mid,eps/2,fl,dep-1)+asr(mid,r,eps/2,fr,dep-1);
}
double calc(double l,double r,double eps){
return asr(l,r,eps,simpson(l,r),12);
}
int main(){
scanf("%lf",&a);
if(a<0){
puts("orz");return 0;
}
double ans=calc(1e-8,20,1e-8);
printf("%.5lf\n",ans);
return 0;
}