C语言实现对称矩阵特征值和特征向量求解

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

#define N 10 // 矩阵最大维数

// 判断矩阵是否为对称矩阵 int isSymmetric(double a[][N], int n) { int i, j; for (i = 0; i < n; i++) { for (j = i + 1; j < n; j++) { if (a[i][j] != a[j][i]) { return 0; } } } return 1; }

// 求矩阵特征值和特征向量 void eigen(double a[][N], int n, double eigval[], double eigvec[][N]) { int i, j, k, m; double b[N][N], c[N][N], p, q, r, s, t, u; double eps = 1e-10; // 精度

// 初始化特征向量矩阵
for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        if (i == j) {
            eigvec[i][j] = 1.0;
        } else {
            eigvec[i][j] = 0.0;
        }
    }
}

// 初始化特征值
for (i = 0; i < n; i++) {
    eigval[i] = a[i][i];
}

// 迭代求解特征值和特征向量
for (m = 0; m < 100; m++) {
    // 判断是否收敛
    double sum = 0.0;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            if (i != j) {
                sum += fabs(a[i][j]);
            }
        }
    }
    if (sum < eps) {
        break;
    }

    // 找到绝对值最大的非对角线元素
    p = q = r = s = 0.0;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            if (i != j) {
                if (fabs(a[i][j]) > fabs(p)) {
                    p = a[i][j];
                    r = i;
                    s = j;
                }
            }
        }
    }

    // 计算旋转角度
    t = (eigval[s] - eigval[r]) / (2.0 * p);
    u = sqrt(1.0 + t * t);
    if (t > 0.0) {
        u = -u;
    }

    // 计算旋转矩阵元素
    b[r][r] = b[s][s] = cos(0.5 * atan(t));
    c[r][s] = -sin(0.5 * atan(t));
    c[s][r] = sin(0.5 * atan(t));
    b[r][s] = -c[r][s];
    for (i = 0; i < n; i++) {
        if (i != r && i != s) {
            b[i][i] = 1.0;
        }
    }
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            eigvec[i][j] = eigvec[i][j] * b[i][j];
        }
    }

    // 计算特征值和特征向量
    eigval[r] = eigval[r] - p;
    eigval[s] = eigval[s] + p;
    for (k = 0; k < n; k++) {
        if (k != r && k != s) {
            p = a[r][k];
            q = a[s][k];
            a[r][k] = b[r][r] * p + b[r][s] * q;
            a[s][k] = b[s][r] * p + b[s][s] * q;
            a[k][r] = a[r][k];
            a[k][s] = a[s][k];
        }
    }
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            if (i != r && i != s && j != r && j != s) {
                eigvec[i][j] = eigvec[i][j];
            } else {
                p = eigvec[i][r];
                q = eigvec[j][s];
                eigvec[i][r] = b[r][r] * p + b[r][s] * q;
                eigvec[j][s] = b[s][r] * p + b[s][s] * q;
            }
        }
    }
}

}

int main() { int n, i, j; double a[N][N], eigval[N], eigvec[N][N];

printf('请输入矩阵维数（最大%d）：', N);
scanf('%d', &n);

printf('请输入矩阵元素：\n');
for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        scanf('%lf', &a[i][j]);
    }
}

if (!isSymmetric(a, n)) {
    printf('输入矩阵不是对称矩阵\n');
    return 0;
}

eigen(a, n, eigval, eigvec);

printf('特征值：\n');
for (i = 0; i < n; i++) {
    printf('%8.4f\n', eigval[i]);
}

printf('特征向量：\n');
for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        printf('%8.4f ', eigvec[i][j]);
    }
    printf('\n');
}

return 0;

}