研究生數值分析上機作業：SOR、列主元消去、三次樣條插值

結果與討論

• 一次編譯通過，關鍵注意：兩重求和變量每輪必歸零；下標 1-based 需統一，避免越界。
• 迭代次數固定 9 次，實際應引入殘差範數 ||Ax-b|| < ε 作為停止條件。

題目一（三）：列主元高斯消去法

理論依據

為提高數值穩定性，每列選擇絕對值最大元素作主元，交換行，消去為上三角矩陣，再回代求解。

核心代碼片段

double a[N+1][N+2]; // 增廣矩陣 [A|b]
// 列主元選取與行交換
for (int j=1; j<=N-1; ++j) {
    int p = j;
    for (int i=j+1; i<=N; ++i)
        if (fabs(a[i][j]) > fabs(a[p][j])) p = i;
    if (p != j) swap_rows(a, j, p, j, N+1); // 從第 j 列開始交換
    // 消去
    for (int i=j+1; i<=N; ++i) {
        double factor = a[i][j] / a[j][j];
        for (int k=j+1; k<=N+1; ++k) // 關鍵：k 從 j+1 開始！
            a[i][k] -= factor * a[j][k];
    }
}
// 回代
x[N] = a[N][N+1] / a[N][N];
for (int i=N-1; i>=1; --i) {
    double sum = 0;
    for (int j=i+1; j<=N; ++j) sum += a[i][j] * x[j];
    x[i] = (a[i][N+1] - sum) / a[i][i];
}

結果與討論

• 致命 Bug：消去內層循環初寫 k=j 導致主元列歸零，後續計算全錯。修正為 k=j+1 即正確。
• 將 A、b 合併為增廣矩陣編程更便捷。

題目三：三次樣條插值 (自然樣條 + 追趕法)

理論依據

區間 [x_i, x_{i+1}] 上構造三次多項式 S_i(x)，滿足：函數值連續、一二階導數連續、邊界二階導數為 0（自然樣條）。導出三對角線性方程組求 M_i = S''(x_i)，再由 M_i 計算插值點函數值與導數。

追趕法求解三對角系統

// 系統：λ_i M_{i-1} + 2 M_i + μ_i M_{i+1} = d_i
// 遞推：P_i = 2 - λ_i * L_{i-1},  L_i = μ_i / P_i,  G_i = (d_i - λ_i * G_{i-1}) / P_i
// 回代：M_n = G_n; M_i = G_i - L_i * M_{i+1}

核心代碼片段

double h = 0.5; // 等距節點
double f[11] = {0, 0, log(2), log(3), ..., log(10)}; // 已知函數值
double lambda[11], mu[11], d[11], M[11];
// 設置 λ, μ, d（邊界條件 M_1=M_10=0 對應 λ_1=μ_10=0）
// 追趕法遞推
double P[11], L[11], G[11];
P[1] = 2; G[1] = d[1];
for (int i=2; i<=9; ++i) {
    L[i] = lambda[i] / P[i-1];
    P[i] = 2 - lambda[i] * L[i];
    G[i] = d[i] - lambda[i] * G[i-1];
}
// 回代
M[10] = G[10] / P[10];
for (int i=9; i>=1; --i) M[i] = (G[i] - mu[i]*M[i+1]) / P[i];
// 計算 x=4.563 處函數值與導數（落在 [4,5] 區間）
double s = M[4]*pow(5-4.563,3)/6 + M[5]*pow(4.563-4,3)/6
         + (f[4]-M[4]/6)*(5-4.563) + (f[5]-M[5]/6)*(4.563-4);
double z = -M[4]/2*pow(5-4.563,2) + M[5]/2*pow(4.563-4,2)
         - (f[4]-M[4]/6) + (f[5]-M[5]/6);

結果與討論

• 下標 0 棄用、1 起始，統一風格可避免混淆。
• 常見錯誤：回代循環 i-- 寫成 i++ 導致死循環；乘法漏寫 *（如 M[4]/2*(...)）。
• 建議引入 assert 檢查分母非零、節點嚴格遞增等前置條件。

總結

三大經典算法完整實現，暴露出下標管理、循環邊界、數值穩定性等工程細節。數值計算「理論正確 ≠ 代碼正確」，嚴謹的邊界檢查、殘差驗證、單元測試不可或缺。