题解 P4525 【【模板】自适应辛普森法1】

· · 题解

写在前面

本题解包含两部分：

  1. 一个我看起来比较舒服的数学推导
  2. 自适应辛普森的一个简单并且没什么卵用的优化

    数学推导

    看好多大佬都用数学方法做的，但是我这个垃圾并看不懂。所以这里给出一个好懂一点的优化。

记被积函数为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

于是就可以积出题目中那个式子了。代码如下：

#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;
}

自适应辛普森

普通的自适应辛普森我就不多说了，代码如下：

#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;
}

然后我翻了一下题解，竟然所有大佬都是这么写的。

仔细看一下这份代码，我们发现，有很多函数值被重复计算了很多次。在本题中，这是不要紧的。但在一些题中，函数值不能用常数时间计算。所以，我们需要维护算过的函数值。代码如下：

#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;
}

我知道这种码风很难看，但是我找不到好看的了……

完结撒花！顺便骗赞