求解特征方程C++

特征方程是指线性代数中矩阵的特征值方程。假设 $A$ 是一个 $n\times n$ 的矩阵，$\lambda$ 是其特征值，那么特征方程为：

$$
\det(A-\lambda I) = 0
$$

其中 $I$ 是 $n\times n$ 的单位矩阵。解特征方程就是求解上述方程的根 $\lambda$。

以下是 C++ 代码实现：

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

const double eps = 1e-10; // 精度

// 求解一元二次方程的根
void solve(double a, double b, double c, vector<double>& res) {
    double delta = b * b - 4 * a * c;
    if (delta < -eps) return; // 无实数解
    if (delta < eps) { // 有一个实根
        double x = -b / (2 * a);
        if (abs(x) < eps) x = 0; // 防止误差
        res.push_back(x);
    } else { // 有两个实根
        double x1 = (-b - sqrt(delta)) / (2 * a);
        if (abs(x1) < eps) x1 = 0; // 防止误差
        res.push_back(x1);
        double x2 = (-b + sqrt(delta)) / (2 * a);
        if (abs(x2) < eps) x2 = 0; // 防止误差
        res.push_back(x2);
    }
}

// 求解特征方程的根
vector<double> characteristic_equation(vector<vector<double>>& A) {
    int n = A.size();
    vector<double> res;
    if (n == 1) { // 1x1 矩阵的特征值就是其唯一元素
        res.push_back(A[0][0]);
        return res;
    }
    vector<double> poly(n + 1); // 保存特征多项式的系数
    poly[n] = 1;
    for (int k = 0; k < n; k++) { // 求解特征多项式
        vector<vector<double>> Ak(n - k, vector<double>(n - k));
        for (int i = k; i < n; i++) {
            for (int j = k; j < n; j++) {
                Ak[i - k][j - k] = A[i][j];
            }
        }
        double a = 1, b = -Ak[0][0];
        for (int i = 1; i < n - k; i++) {
            a = -Ak[i][i] * a;
            b = -Ak[i][i] * b - Ak[i - 1][i] * a;
        }
        solve(a, b, -Ak[n - k - 1][n - k - 1], res);
        poly[n - k - 1] = -b / a;
    }
    return res;
}

int main() {
    vector<vector<double>> A = {{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}};
    vector<double> res = characteristic_equation(A);
    for (auto x : res) {
        cout << x << " ";
    }
    cout << endl;
    return 0;
}

上述代码实现了求解特征方程和特征值的功能。