一、复化辛卜生求积公式推导

1. 辛卜生公式基础

辛卜生公式是基于二次插值多项式的求积公式。对于单个区间 [a, b]，取中点 ，构造过 (a, f(a))、(c, f(c))、(b, f(b)) 的二次插值多项式 ，则：
通过计算二次插值积分，得到单区间辛卜生公式： 其中 h = b - a（区间长度）。

2. 复化辛卜生公式构造

为提高精度，将积分区间 [a, b] n 等分（n 必须为偶数），步长 ，分割点为：
在每个子区间 （长度 2h）上应用辛卜生公式，总积分近似为所有子区间积分之和：
整理后得到复化辛卜生公式：

• 第一项 f(a) 和最后一项 f(b)：端点函数值，系数为 1；
• 中间奇数项 f(x_{2k+1})：子区间中点函数值，系数为 4；
• 中间偶数项 f(x_{2k})：子区间分点函数值，系数为 2。

二、变步长复化辛卜生求积算法

#include <stdio.h>  
#include <math.h>  
  
double simpson(double (*f)(double x), double a, double b, double e, int max) {  
    int i, k;  
    double h = b - a;  
    // 初始梯形积分 T0 (将区间分为1个子区间)  
    double T_old = h * (f(a) + f(b)) / 2.0;  
    double T_new, S_old, S_new;  
  
    // 第一次分半：计算T1（2个子区间的梯形积分）和S1（2个子区间的辛普森积分）  
    h /= 2.0; // 步长减半  
    T_new = T_old / 2.0 + h * f(a + h); // 新增1个中点  
    S_old = (4.0 * T_new - T_old) / 3.0; // 辛普森积分公式  
  
    if (max <= 1) return S_old;  
  
    // 逐次分半迭代  
    for (i = 1; i < max; i++) {  
        h /= 2.0; // 步长再次减半  
        double sum_mid = 0.0;  
        // 计算本次新增的所有中点（数量=2^i）  
        for (k = 0; k < (1 << i); k++) { // 1<<i 等价于 2^i            double x = a + (2 * k + 1) * h; // 新增中点的坐标  
            sum_mid += f(x);  
        }  
        // 递推更新梯形积分  
        T_old = T_new;  
        T_new = T_old / 2.0 + h * sum_mid;  
        // 由梯形积分计算新的辛普森积分  
        S_new = (4.0 * T_new - T_old) / 3.0;  
  
        // 误差满足则返回结果  
        if (fabs(S_new - S_old) < e) {  
            return S_new;  
        }  
  
        S_old = S_new;  
    }  
  
    return S_old;  
}  

三、验证实例

1. 验证用例选择

选择已知解析解的积分进行验证，便于对比结果：

解析解：

2. 完整验证代码

#include <stdio.h>  
#include <math.h>  
  
double simpson(double (*f)(double x), double a, double b, double e, int max) {  
    int i, k;  
    double h = b - a;  
    // 初始梯形积分 T0 (将区间分为1个子区间)  
    double T_old = h * (f(a) + f(b)) / 2.0;  
    double T_new, S_old, S_new;  
  
    // 第一次分半：计算T1（2个子区间的梯形积分）和S1（2个子区间的辛普森积分）  
    h /= 2.0; // 步长减半  
    T_new = T_old / 2.0 + h * f(a + h); // 新增1个中点  
    S_old = (4.0 * T_new - T_old) / 3.0; // 辛普森积分公式  
  
    if (max <= 1) return S_old;  
  
    // 逐次分半迭代  
    for (i = 1; i < max; i++) {  
        h /= 2.0; // 步长再次减半  
        double sum_mid = 0.0;  
        // 计算本次新增的所有中点（数量=2^i）  
        for (k = 0; k < (1 << i); k++) { // 1<<i 等价于 2^i            double x = a + (2 * k + 1) * h; // 新增中点的坐标  
            sum_mid += f(x);  
        }  
        // 递推更新梯形积分  
        T_old = T_new;  
        T_new = T_old / 2.0 + h * sum_mid;  
        // 由梯形积分计算新的辛普森积分  
        S_new = (4.0 * T_new - T_old) / 3.0;  
  
        // 误差满足则返回结果  
        if (fabs(S_new - S_old) < e) {  
            return S_new;  
        }  
  
        S_old = S_new;  
    }  
  
    return S_old;  
}  
  
// 测试函数：f(x) = 1/(1+x³)  
double test_func(double x) {  
    return 1.0 / (1.0 + x * x * x);  
}  
  
int main() {  
    double a = 0.0, b = 1.0;  
    double epsilon = 1e-8;
    int max_iter = 20;  
  
    double result = simpson(test_func, a, b, epsilon, max_iter);  
    printf("积分结果: %.8f\n", result); // 输出：0.83564885  
  
    return 0;  
}

3. 运行结果与分析

• 输出结果：变步长复化辛卜生求积结果：0.83564885
• 误差分析：计算结果与解析解的误差小于 10^{-6}，满足精度要求；