对于一个二次函数 \(f(x) = ax^2 + bx + c\)，积分得 \(F(x) = \displaystyle\int_0^x f(t) \, \mathrm{d}t = \dfrac{a}{3}x^3 + \dfrac{b}{2}x^2 + cx + C\)。
于是
\[\displaystyle\int_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} \]以上即为辛普森公式。
普通辛普森法
为了求 \(\displaystyle\int^r_l f(x) \, \mathrm{d}x\)，可以将区间 \([l, r]\) 分为\(2n\)个等长的区间 \([g_i, g_{i + 1}](i \in \{1, 2, ..., 2n\})\)。设 \(h = \dfrac{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^r_l f(x) \, \mathrm{d}x \approx \displaystyle\sum_{i = 1}^{n} \displaystyle\int^{g_{2i + 1}}_{g_{2i - 1}} P(x) \, \mathrm{d}x \]\[= \dfrac{h}{3} (\displaystyle\sum^{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;
}
自适应辛普森法
考虑到当分割出的二次函数在该区间内已经足够接近原函数时，可以直接使用该积分值；而当二次函数仍不够接近时还需再分割，于是可以考虑分别用辛普森公式计算区间 \([l, r]\)、\([l, \dfrac{l + r}{2}]\)、\([\dfrac{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;
}