C语言Simpson积分怎么实现？

辛普森积分法的原理

辛普森法基于一个数学结论：对于任意一个不超过三次的多项式函数 f(x)，在区间 [a, b] 上的定积分，等于该区间中点 c = (a+b)/2 处的函数值 f(c) 加上两个端点处的函数值 f(a) 和 f(b)，再乘以一个系数 (b-a)/6。

即： $$ \int_a^b f(x) \,dx \approx \frac{b-a}{6} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right] $$

这个公式被称为辛普森 1/3 规则。

为了提高精度,我们可以将整个积分区间 [a, b] 分成 n 个（n 必须是偶数）等宽的小子区间，然后在每个小子区间上应用上述 1/3 规则，最后将所有结果相加。

如果我们将区间 [a, b] 分成 n 个等宽的子区间，每个子区间的宽度为 h = (b - a) / n，那么积分可以近似为：

$$ \int_a^b f(x) \,dx \approx \frac{h}{3} \left[ f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x3) + \dots + 2f(x{n-2}) + 4f(x_{n-1}) + f(x_n) \right] $$

x_i = a + i * h。

系数模式：1, 4, 2, 4, 2, ..., 4, 1，这个模式是辛普森法的核心。

C 语言实现

下面我们分步来实现这个算法。

步骤 1：定义要积分的函数

我们需要一个函数来表示我们要计算积分的数学函数 f(x)，我们可以用一个 C 函数指针来指向它，这样我们的积分函数就可以处理任何符合该签名的函数。

// 示例函数1: f(x) = x^2
double f1(double x) {
    return x * x;
}
// 示例函数2: f(x) = sin(x)
#include <math.h> // 需要包含 math.h 来使用 sin()
double f2(double x) {
    return sin(x);
}

步骤 2：编写辛普森积分函数

这是核心部分,我们将创建一个函数 simpson_integration，它接收函数指针、积分下限 a、积分上限 b 和子区间数量 n。

#include <stdio.h>
#include <math.h>
// 函数指针类型，指向一个接受 double 返回 double 的函数
typedef double (*FuncPtr)(double);
/**
 * @brief 使用辛普森1/3规则计算定积分
 * 
 * @param func 指向被积函数的指针
 * @param a 积分下限
 * @param b 积分上限
 * @param n 子区间的数量 (必须是偶数)
 * @return double 积分的近似值
 */
double simpson_integration(FuncPtr func, double a, double b, int n) {
    // 1. 检查 n 是否为偶数
    if (n % 2 != 0) {
        printf("错误: 子区间数量 n 必须是偶数，\n");
        return 0.0; // 或者可以返回一个错误码
    }
    // 2. 计算步长 h
    double h = (b - a) / n;
    // 3. 初始化求和结果
    double integral = func(a) + func(b); // f(x0) + f(xn)
    // 4. 循环计算中间项
    for (int i = 1; i < n; i++) {
        double x = a + i * h;
        if (i % 2 == 1) {
            // 奇数索引项: 4 * f(x_i)
            integral += 4.0 * func(x);
        } else {
            // 偶数索引项: 2 * f(x_i)
            integral += 2.0 * func(x);
        }
    }
    // 5. 乘以 h/3 得到最终结果
    integral = integral * h / 3.0;
    return integral;
}

步骤 3：主函数 main - 完整示例

我们把所有部分组合起来,在 main 函数中进行测试。

#include <stdio.h>
#include <math.h>
// 函数指针类型
typedef double (*FuncPtr)(double);
// 定义要积分的函数
double f1(double x) {
    return x * x;
}
double f2(double x) {
    return sin(x);
}
// 辛普森积分函数
double simpson_integration(FuncPtr func, double a, double b, int n) {
    if (n % 2 != 0) {
        printf("错误: 子区间数量 n 必须是偶数，\n");
        return 0.0;
    }
    double h = (b - a) / n;
    double integral = func(a) + func(b);
    for (int i = 1; i < n; i++) {
        double x = a + i * h;
        if (i % 2 == 1) {
            integral += 4.0 * func(x);
        } else {
            integral += 2.0 * func(x);
        }
    }
    return integral * h / 3.0;
}
int main() {
    // --- 测试 1: 计算 ∫(0 to 2) x^2 dx ---
    // 精确解: [x^3/3] from 0 to 2 = 8/3 ≈ 2.666666...
    double a1 = 0.0, b1 = 2.0;
    int n1 = 100; // 使用100个子区间
    double result1 = simpson_integration(f1, a1, b1, n1);
    printf("计算 ∫(%.1f 到 %.1f) x^2 dx:\n", a1, b1);
    printf("辛普森法结果 (n=%d): %.8f\n", n1, result1);
    printf("精确解: 8/3 ≈ %.8f\n\n", 8.0 / 3.0);
    // --- 测试 2: 计算 ∫(0 to π) sin(x) dx ---
    // 精确解: [-cos(x)] from 0 to π = -cos(π) - (-cos(0)) = -(-1) - (-1) = 2
    double a2 = 0.0, b2 = M_PI; // M_PI 是 math.h 中定义的 π
    int n2 = 100;
    double result2 = simpson_integration(f2, a2, b2, n2);
    printf("计算 ∫(%.1f 到 %.1f) sin(x) dx:\n", a2, b2);
    printf("辛普森法结果 (n=%d): %.8f\n", n2, result2);
    printf("精确解: 2\n");
    return 0;
}

编译和运行

你需要使用 C 编译器（如 GCC）来编译这个程序，由于我们使用了 math.h 中的 sin 和 M_PI，需要链接数学库。

在 Linux/macOS 上编译：

gcc simpson.c -o simpson -lm

-lm 选项告诉链接器链接数学库。

在 Windows (MinGW) 上编译：

gcc simpson.c -o simpson.exe -lm

运行程序：

./simpson

预期输出：

计算 ∫(0.0 到 2.0) x^2 dx:
辛普森法结果 (n=100): 2.66666667
精确解: 8/3 ≈ 2.66666667
计算 ∫(0.0 到 3.141593) sin(x) dx:
辛普森法结果 (n=100): 2.00000000
精确解: 2

可以看到,辛普森法即使只用了 100 个子区间，也能得到非常精确的结果。

代码优化与变体

上面的实现是清晰易懂的,但在某些情况下可以进行优化。

优化：减少 func() 调用次数

上面的实现在循环中调用了 n-1 次 func()，我们可以通过重用计算结果来减少函数调用次数，但这会增加代码的复杂性，对于大多数应用场景，可读性比微小的性能提升更重要，因此原版是更好的选择。

变体：复合辛普森 3/8 规则

当子区间数量 n 是3的倍数时，可以使用辛普森 3/8 规则，它在每个由三个子区间组成的区间上使用三次多项式来近似。

$$ \int_a^b f(x) \,dx \approx \frac{3h}{8} \left[ f(x_0) + 3f(x_1) + 3f(x_2) + 2f(x3) + \dots + 3f(x{n-2}) + 3f(x_{n-1}) + f(x_n) \right] $$

h = (b-a)/n。

系数模式：1, 3, 3, 2, 3, 3, 2, ..., 3, 3, 1。

你可以基于这个规则编写一个类似的 simpson_3_8_integration 函数。

• 核心思想：用抛物线近似函数曲线，比梯形法更精确。
• 关键公式：integral = (h/3) * (f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + ... + f(xn))
• C 语言实现要点：
  1. 使用函数指针 (FuncPtr) 使积分函数通用化。
  2. 确保 n 是偶数。
  3. 循环中根据索引的奇偶性应用不同的系数（4 或 2）。
  4. 编译时记得链接数学库 (-lm)。
• 优点：精度高，对于平滑函数收敛速度很快。
• 缺点：要求 n 为偶数，对于某些函数（如有剧烈振荡或奇点的函数），可能需要非常多的 n 才能收敛。