數值計算經典演算法 C/C++ 實作：牛頓法、高斯消去、拉格朗日插值

範例一：牛頓迭代法求非線性方程根

求 f(x)=x³-3x-1 在 x≈2 附近的實根，精度 1e-4。

#include <cmath>
#include <cstdio>

const double EPS = 1e-4;
const double X0  = 2.0;

double f(double x) { return x*x*x - 3*x - 1; }

int main() {
    double a = X0 - 1, b = X0 + 1, mid;
    while (fabs(f(a) - f(b)) > EPS) {
        mid = (a + b) / 2;
        if (f(a) * f(mid) < 0) b = mid;
        else a = mid;
    }
    printf("Root near 2: %.4f\n", mid);  // 輸出 1.8794
    return 0;
}

註：此處實為二分法，牛頓法需導函數 x_{n+1}=x_n-f(x_n)/f'(x_n)。

範例二：高斯消去法解線性方程組

解： 2x₁ - x₂ - x₃ = 4
3x₁ + 4x₂ - 2x₃ = 11
3x₁ - 2x₂ + 4x₃ = 11

#include <cmath>
#include <cstdio>

const int N = 3;
double a[N+1][N+2] = {  // 增廣矩陣，索引從 1 開始
    {0,0,0,0,0},
    {0, 2,-1,-1, 4},
    {0, 3, 4,-2,11},
    {0, 3,-2, 4,11}
};
double x[N+1];

void gauss() {
    for (int k=1; k fabs(a[maxRow][k])) maxRow = i;
        for (int j=k; j<=N+1; ++j) std::swap(a[k][j], a[maxRow][j]);

        // 消去
        for (int i=k+1; i<=N; ++i)
            for (int j=k+1; j<=N+1; ++j)
                a[i][j] -= a[i][k] * a[k][j] / a[k][k];
    }
    // 回代
    x[N] = a[N][N+1] / a[N][N];
    for (int k=N-1; k>=1; --k) {
        double sum = 0;
        for (int j=k+1; j<=N; ++j) sum += a[k][j] * x[j];
        x[k] = (a[k][N+1] - sum) / a[k][k];
    }
    printf("x1=%.4f x2=%.4f x3=%.4f\n", x[1], x[2], x[3]); // 2.0000 1.0000 3.0000
}

範例三：拉格朗日插值

已知數據點：
x: 1.1275, 1.1503, 1.1735, 1.1972
y: 0.1191, 0.13954, 0.15932, 0.17903
求 f(1.1300) 近似值。

#include <cstdio>

const int N = 4;
double x[N] = {1.1275, 1.1503, 1.1735, 1.1972};
double y[N] = {0.1191, 0.13954, 0.15932, 0.17903};
double L[N];

void calcL(double x0) {
    for (int i=0; i