C 語言數值積分實作：六大經典演算法完整解析

前言

定積分數值求解是科學計算的基石。本文整理六種常見演算法——蒙地卡羅隨機點法、左矩形法、梯形法、復合梯形法、辛普生法、高斯求積法——並提供完整、可編譯的 C 語言實作碼，適合教學參考與工程移植。

核心函數原型

double fun(double x);           // 被積函數
double Darts(int n);            // 隨機點法
double LeftRect(double a, double b, int n);  // 左矩形法
double Trape(double a, double b, int n);     // 梯形法
double T(double a, double b, int n);         // 復合梯形法
double Simposn(double a, double b, int n);   // 辛普生法
double Gass(double (*f)(double), double a, double b, int n); // 高斯求積法

演算法細節與程式碼片段

1. 蒙地卡羅隨機點法

利用亂數在單位正方形內投點，統計落在曲線下的比例估算積分值：

double Darts(int n) {
    double x, y; int count = 0;
    srand((unsigned)time(NULL));
    for (int i = 0; i < n; i++) {
        x = rand() % 100 / 100.0;
        y = rand() % 100 / 100.0;
        if (y <= 1 - x * x) count++;
    }
    return (double)count / n;
}

2. 左矩形法

區間等分，以各子區間左端點函數值為高構成矩形求和：

double LeftRect(double down, double up, int n) {
    double h = (up - down) / n;
    double sum = fun(down) * h;
    for (int i = 1; i < n; i++)
        sum += fun(down + i * h) * h;
    return sum;
}

3. 梯形法與復合梯形法

梯形法以兩端點連線近似；復合梯形法進一步細分區間，精度更高：

double Trape(double down, double up, int n) {
    double h = (up - down) / n;
    double sum = 0.5 * (fun(down) + fun(down + h)) * h;
    for (int i = 1; i < n; i++)
        sum += 0.5 * (fun(down + i * h) + fun(down + (i + 1) * h)) * h;
    return sum;
}

double T(double x, double y, int z) {  // 復合梯形
    double h = (y - x) / z;
    double Tn = (fun(x) + fun(y)) / 2.0;
    for (int i = 1; i < z; i++) Tn += fun(x + i * h);
    return Tn * h;
}

4. 辛普生法

利用拋物線插值，公式 S = (4*T(2n) - T(n)) / 3，精度顯著優於梯形法：

double Simposn(double down, double up, int n) {
    return (4 * T(down, up, 2 * n) - T(down, up, n)) / 3.0;
}

5. 高斯求積法

預設 5 階高斯點與權重，針對光滑函數極高效：

double Gass(double (*func)(double), double a, double b, int n) {
    static const double *nodes[] = { /* x1..x5 */ };
    static const double *weights[] = { /* A1..A5 */ };
    double sum = 0;
    for (int i = 0; i < n; i++)
        sum += func((b - a) * nodes[n-1][i] / 2 + (a + b) / 2) * weights[n-1][i];
    return sum * (b - a) / 2;
}

主程式測試流程

程式依序呼叫各演算法，並輸出結果供比對：

int main() {
    printf("隨機點法: %f\n", Darts(10000));
    double a, b; int n;
    scanf("%lf %lf %d", &a, &b, &n);
    printf("左矩形: %f\n", LeftRect(a, b, n));
    printf("梯形法: %f\n", Trape(a, b, n));
    printf("辛普生: %f\n", Simposn(a, b, n));
    printf("高斯法: %f\n", Gass(fun, a, b, 3));
    return 0;
}

效能與精度比較建議

常見坑點提醒

• rand() % 100 / 100.0 僅產生 0.00~0.99，若需 [0,1) 應改用 rand() / (RAND_MAX + 1.0)。
• 辛普生法要求 n 為偶數，否則公式不成立。
• 高斯求積的節點與權重表僅內建至 5 階，更高階需自行擴充或引用數學庫。