辛普森公式

这里用一个抛物线求弧长的例子给出Simpson公式的代码：

#include <cstdio>
#include <cmath>
const double eps = 1e-5;
// 被积函数F
double a;
double F(double x){
    return sqrt(1 + 4 * a * a * x * x);
}
// 三点Simpson法。这里要求F是一个全局函数，F也就是被积函数
double simpson(double a,double b){
    double c = a + (b - a) / 2;
    return (F(a) + 4 * F(c) + F(b)) * (b - a) / 6;
}
// 自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点Simpson值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.0;
    return asr(a,c,eps / 2,L) + asr(c,b,eps / 2,R);
}
// 自适应Simpson公式（主过程）
double asr(double a,double b,double eps){
    return asr(a,b,eps,simpson(a,b));
}
// 用自适应Simpson公式计算宽为w
double parabola_arc_length(double w,double h){
    a = 4.0 * h / (w * w); // 修改全局变量a，从而影响F的行为
    return asr(0,w / 2,eps) * 2;
}