C語言實現預估校正法求解常微分方程

C語言實現預估校正法求解常微分方程

預估校正法是求解常微分方程初值問題的經典數值方法之一。本文將通過具體的C語言代碼示例，詳細講解如何使用預估校正法（Predictor-Corrector Method）來求解常微分方程 y' = 1 + y²，初始條件為 y(0) = 0，在區間 [0, 1.4] 上的數值解。

算法原理

預估校正法屬於多步法，其核心思想分為兩步：

• 預估步驟：使用顯式公式（如歐拉法）預估下一步的近似值
• 校正步驟：使用隱式公式（如梯形公式）利用預估值進行修正，獲得更高精度的結果

本例採用歐拉法作為預估器，梯形公式作為校正器，具體公式為：

預估：y_{i+1}^{(0)} = y_i + h * f(x_i, y_i)
校正：y_{i+1} = y_i + h/2 * [f(x_i, y_i) + f(x_{i+1}, y_{i+1}^{(0)})]

完整C語言實現代碼

#include <stdio.h>
#include <math.h>

#define N 14

int main() {
    double a = 0.0, b = 1.4, h, m = 0.0;
    double x[N + 1] = {0.0}, y[N + 1] = {0};

    h = (b - a) / N;
    x[0] = a;
    y[0] = m;

    for (int i = 1; i <= N; i++) {
        // 預估步驟：歐拉法
        y[i] = y[i - 1] + h * (1 + y[i - 1] * y[i - 1]);
        x[i] = a + i * h;
        // 校正步驟：梯形公式
        y[i] = y[i - 1] + h * 0.5 * ((1 + y[i - 1] * y[i - 1]) + (1 + y[i] * y[i]));
        printf("x[%d]=%e , y[%d]=%e\n", i, x[i], i, y[i]);
    }
    return 0;
}

代碼關鍵點解析

• 步長計算：h = (b - a) / N 將區間等分為 N 段
• 數組存儲：使用 x[] 和 y[] 分別存儲自變量和因變量的離散值
• 疊代過程：每次循環先用歐拉法得到預估值，再用梯形公式校正
• 微分方程：f(x, y) = 1 + y² 直接嵌入在計算公式中

運行結果示例

程序輸出了各網格點的 x 和 y 值，展示了數值解隨自變量變化的趨勢。由於該微分方程解析解為 y = tan(x)，在 x 接近 π/2 ≈ 1.57 時發散，因此在 1.4 處數值仍保持穩定。

總結

預估校正法結合了顯式方法的簡單性和隱式方法的穩定性，是工程計算中常用的數值算法。通過調整步長 N 和選擇不同的預估/校正公式組合，可以在精度和計算量之間取得平衡。