自适应辛普森法 学习笔记

对于一个二次函数 $f(x)=a{x}^2+bx+c$，积分得 $F(x)={\displaystyle {\int }_0^{x}f(t)\mathrm{d}t={\displaystyle \frac{a}{3}}{x}^3+{\displaystyle \frac{b}{2}}{x}^2+cx+C}$。

于是

$${\displaystyle {\int }_l^{r}f(x)\mathrm{d}x=F(r)-F(l)}$$

$$={\displaystyle \frac{a(r^3-l^3)}{3}}+{\displaystyle \frac{b(r^2-l^2)}{2}}+c(r-l)+(C-C)$$

$$=(r-l)[{\displaystyle \frac{a(r^2+rl+l^2)}{3}}+{\displaystyle \frac{b(r+l)}{2}}+c]$$

$$={\displaystyle \frac{r-l}{6}}(2a{r}^2+2arl+2a{l}^2+3br+3bl+6c)$$

$$={\displaystyle \frac{r-l}{6}}\{(a{r}^2+br+c)+(a{l}^2+bl+c)+4[a({\displaystyle \frac{l+r}{2}}{)}^2+b({\displaystyle \frac{l+r}{2}})+c]\}$$

$$={\displaystyle \frac{(r-l)(f(l)+4f({\displaystyle \frac{l+r}{2}})+f(r))}{6}}$$

以上即为辛普森公式。

普通辛普森法

为了求 ${\displaystyle {\int }_l^{r}f(x)\mathrm{d}x}$，可以将区间 $[l,r]$ 分为$2n$个等长的区间 $[g_i,g_{i+1}](i\in \{1,2,...,2n\})$。设 $h={\displaystyle \frac{r-l}{2n}}$，近似地计算过 $(g_{i-1},f(g_{i-1}))$、$(g_i,f(g_i))$、$(g_{i+1},f(g_{i+1}))$ 三点的二次函数在区间 $[g_{i-1},g_{i+1}]$ 的积分。当 $n$ 足够大时，计算结果便越接近原函数的积分。设 $P(x)$ 为过该三点的函数，则

$${\displaystyle {\int }_l^{r}f(x)\mathrm{d}x\approx {\displaystyle \sum _{i=1}^n{\displaystyle {\int }_{g_{2i-1}}^{g_{2i+1}}P(x)\mathrm{d}x}}}$$

$$={\displaystyle \frac{h}{3}}({\displaystyle \sum _{i=1}^{n}f(g_{2i-1})+4f(g_{2i})+f(g_{2i+1}))}$$

$$={\displaystyle \frac{h}{3}}(f(g_1)+4f(g_2)+2f(g_3)+4f(g_4)+...+2f(g_{2n-1})+4f(g_{2n})+f(g_{2n+1}))$$

Code:

int n = 1e7;

double calc(double l, double r) {
    double h = (r - l) / (2 * n), res = f(l) + f(r);
    
    for (int i = 2; i < 2 * n; ++i)
        res += f(l + i * h) * (i & 1? 2: 4);
    
    return res * h / 3;
}

自适应辛普森法

考虑到当分割出的二次函数在该区间内已经足够接近原函数时，可以直接使用该积分值；而当二次函数仍不够接近时还需再分割，于是可以考虑分别用辛普森公式计算区间 $[l,r]$、$[l,{\displaystyle \frac{l+r}{2}}]$、$[{\displaystyle \frac{l+r}{2}},r]$ 的积分值，若大区间的积分值与小区间的积分值和的差足够小，则可以认为二次函数的图像在该区间内足够接近原函数的图像，并直接返回大区间的积分值。否则，递归计算两个小区间。一般还需设定一个次数，当已经递归超过该次数时才开始判断是否可以直接返回。

P4525 【模板】自适应辛普森法 1 AC Code:

#include <cstdio>
#define simpson(l, r) ((r - l) * (f(l) + f(r) + 4 * f((l + r) / 2)) / 6)
#define f(x) (double(c * (x) + d) / (a * (x) + b))
#define abs(x) ((x) >= 0? (x): -(x))

using namespace std;

double a, b, c, d, L, R;

double calc(double l, double r, int step, double sps) {
    double mid = (l + r) / 2;
    double fl = simpson(l, mid), fr = simpson(mid, r);
    
    if (abs(fl + fr - sps) <= 1e-6 && step < 0)
        return sps;
    else
        return calc(l, mid, step - 1, fl) + calc(mid, r, step - 1, fr);
}

int main() {
    scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &L, &R);
    printf("%lf", calc(L, R, 12, simpson(L, R)));
    
    return 0;
}