题解 P4525 【【模板】自适应辛普森法1】
FZzzz
2019-07-14 12:34:03
## 写在前面
本题解包含两部分：
1. 一个我看起来比较舒服的数学推导
1. 自适应辛普森的一个简单~~并且没什么卵用~~的优化
## 数学推导
看好多大佬都用数学方法做的，但是我这个垃圾并看不懂。所以这里给出一个好懂一点的优化。
记被积函数为$f(x)$。
首先，我们谈判$a=0$的情况。
这种情况比较简单，我就不说了。（逃
若$a\not=0$，将被积函数进行代数变形
$$\begin{aligned}f(x)&=\frac{cx+d}{ax+b}\\&=\frac{c}{a}+\frac{\frac{ad-bc}{a}}{ax+b}\end{aligned}$$
则
$$\int f(x)\mathrm{d}x=\int\frac{c}{a}\mathrm{d}x+\frac{ad-bc}{a}\int\frac{\mathrm{d}x}{ax+b}$$
加号左边那个积分是很容易积出来的，而查一下积分表，有
$$\int\frac{\mathrm{d}x}{ax+b}=\frac{1}{a}\ln|ax+b|+C$$
故
$$\int f(x)\mathrm{d}x=\frac{ad-bc}{a^2}\ln|ax+b|+\frac{cx}{a}+C$$
于是就可以积出题目中那个式子了。代码如下：
```cpp
#include<cmath>
#include<cstdio>
double a,b,c,d;
double F(double x){
return (a*d-b*c)*log(fabs(a*x+b))/(a*a)+c*x/a;
}
double F2(double x){
return c*x*x/(2*b)+d*x/b;
}
int main(){
#ifdef LOCAL
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
double l,r;
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
if(a!=0.0) printf("%.6f\n",F(r)-F(l));
else printf("%.6f\n",F2(r)-F2(l));
return 0;
}
```
## 自适应辛普森
普通的自适应辛普森我就不多说了，代码如下：
```cpp
#include<cmath>
#include<cstdio>
double a,b,c,d;
double f(double x){
return (c*x+d)/(a*x+b);
}
//区间[a,b]上的辛普森值
double simpson(double a,double b){
double c=a+(b-a)/2;
return (f(a)+4*f(c)+f(b))*(b-a)/6;
}
//区间[a,b]上的积分，精度限制为eps，已知整个区间的辛普森值A
double asr(double a,double b,double eps,double A){
double c=a+(b-a)/2;
double L=simpson(a,c),R=simpson(c,b);
if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15;
else return asr(a,c,eps/2,L)+asr(c,b,eps/2,R);
}
const double eps=1e-7;
int main(){
#ifdef LOCAL
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
double l,r;
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
printf("%.6f\n",asr(l,r,eps,simpson(l,r)));
return 0;
}
```
然后我翻了一下题解，竟然所有大佬都是这么写的。
仔细看一下这份代码，我们发现，有很多函数值被重复计算了很多次。在本题中，这是不要紧的。但在一些题中，函数值不能用常数时间计算。所以，我们需要维护算过的函数值。代码如下：
```cpp
#include<cmath>
#include<cstdio>
double a,b,c,d;
double f(double x){
return (c*x+d)/(a*x+b);
}
//区间[a,b]上的辛普森值，已知区间长度l和端点及中点处的函数值A，B，C
double simpson(double l,double A,double B,double C){
return (A+4*C+B)*l/6;
}
double asr(double a,double b,double eps,double A,double B,double C){
double l=b-a;
double c=a+l/2;
double D=f(a+l/4),E=f(a+3*l/4);
double L=simpson(l/2,A,C,D),R=simpson(l/2,C,B,E),AB=simpson(l,A,B,C);
if(fabs(L+R-AB)<=15*eps) return L+R+(L+R-AB)/15;
else return asr(a,c,eps/2,A,C,D)+asr(c,b,eps/2,C,B,E);
}
const double eps=1e-7;
int main(){
#ifdef LOCAL
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
double l,r;
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
printf("%.6f\n",asr(l,r,eps,f(l),f(r),f(l+(r-l)/2)));
return 0;
}
```
我知道这种码风很难看，但是我找不到好看的了……
------------
完结撒花！~~顺便骗赞~~