• 欢迎光临~

# 自适应辛普森法 学习笔记

[displaystyleint_l^r f(x) , mathrm{d}x = F(r) - F(l) ]

[= dfrac{a(r^3 - l^3)}{3} + dfrac{b(r^2 - l^2)}{2} + c(r - l) + (C - C) ]

[= (r - l)[dfrac{a(r^2 + rl + l^2)}{3} + dfrac{b(r + l)}{2} + c] ]

[= dfrac{r - l}{6}(2ar^2 + 2arl + 2al^2 + 3br + 3bl + 6c) ]

[= dfrac{r - l}{6}{(ar^2 + br + c) + (al^2 + bl + c) + 4[a(dfrac{l + r}{2})^2 + b(dfrac{l + r}{2}) + c]} ]

[= dfrac{(r - l)(f(l) + 4f(dfrac{l + r}{2}) + f(r))}{6} ]

## 普通辛普森法

[displaystyleint^r_l f(x) , mathrm{d}x approx displaystylesum_{i = 1}^{n} displaystyleint^{g_{2i + 1}}_{g_{2i - 1}} P(x) , mathrm{d}x ]

[= dfrac{h}{3} (displaystylesum^{n}_{i = 1} f(g_{2i - 1}) + 4f(g_{2i}) + f(g_{2i + 1})) ]

[= dfrac{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;
}
``````

## 自适应辛普森法

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