求三阶矩阵的特征值与特征向量c语言

好的，以下是一个求解三阶矩阵特征值与特征向量的 C 语言代码示例：

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

#define N 3

void matrix_multiply(double a[N][N], double b[N][N], double c[N][N]) {
    // 矩阵乘法
    int i, j, k;
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            c[i][j] = 0;
            for (k = 0; k < N; k++) {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
}

void matrix_transpose(double a[N][N], double b[N][N]) {
    // 矩阵转置
    int i, j;
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            b[i][j] = a[j][i];
        }
    }
}

void matrix_print(double a[N][N]) {
    // 打印矩阵
    int i, j;
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%f ", a[i][j]);
        }
        printf("\n");
    }
}

void eigenvalues(double a[N][N], double lambda[N]) {
    // 求特征值
    double A[N][N], B[N][N], C[N][N], D[N][N];
    matrix_multiply(a, a, A);
    matrix_multiply(a, A, B);
    matrix_multiply(A, A, C);
    matrix_multiply(A, B, D);

    double a2 = A[0][0] + A[1][1] + A[2][2];
    double a1 = B[0][0] + B[1][1] + B[2][2];
    double a0 = D[0][0] + D[1][1] + D[2][2];

    double p = a1 / 3.0 - a2 * a2 / 9.0;
    double q = a2 * a1 / 6.0 - a0 / 2.0 - a2 * a2 * a2 / 27.0;
    double delta = p * p * p + q * q;

    if (delta > 0) {
        double sqrt_delta = sqrt(delta);
        double u = pow(-q + sqrt_delta, 1.0 / 3.0);
        double v = pow(-q - sqrt_delta, 1.0 / 3.0);
        lambda[0] = u + v - a2 / 3.0;
    } else if (delta == 0) {
        double u = pow(-q, 1.0 / 3.0);
        lambda[0] = 2.0 * u - a2 / 3.0;
        lambda[1] = -u - a2 / 3.0;
    } else {
        double sqrt_delta = sqrt(-delta);
        double rho = pow(sqrt_delta + fabs(q), 1.0 / 3.0);
        double theta = atan(sqrt_delta / q);
        double t = 2.0 * sqrt(-p / 3.0);
        lambda[0] = t * cos(theta / 3.0) - a2 / 3.0;
        lambda[1] = t * cos((theta + 2.0 * M_PI) / 3.0) - a2 / 3.0;
        lambda[2] = t * cos((theta - 2.0 * M_PI) / 3.0) - a2 / 3.0;
    }
}

void eigenvectors(double a[N][N], double lambda[N], double v[N][N]) {
    // 求特征向量
    double A[N][N], B[N][N], C[N][N];
    int i;
    for (i = 0; i < N; i++) {
        A[i][0] = a[i][0] - lambda[0];
        A[i][1] = a[i][1];
        A[i][2] = a[i][2];

        B[i][0] = a[i][0];
        B[i][1] = a[i][1] - lambda[0];
        B[i][2] = a[i][2];

        C[i][0] = a[i][0];
        C[i][1] = a[i][1];
        C[i][2] = a[i][2] - lambda[0];
    }

    double A_det = A[0][0] * (A[1][1] * A[2][2] - A[2][1] * A[1][2])
            - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0])
            + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
    double B_det = B[0][0] * (B[1][1] * B[2][2] - B[2][1] * B[1][2])
            - B[0][1] * (B[1][0] * B[2][2] - B[1][2] * B[2][0])
            + B[0][2] * (B[1][0] * B[2][1] - B[1][1] * B[2][0]);
    double C_det = C[0][0] * (C[1][1] * C[2][2] - C[2][1] * C[1][2])
            - C[0][1] * (C[1][0] * C[2][2] - C[1][2] * C[2][0])
            + C[0][2] * (C[1][0] * C[2][1] - C[1][1] * C[2][0]);

    if (A_det != 0) {
        double A_inv[N][N];
        A_inv[0][0] = (A[1][1] * A[2][2] - A[2][1] * A[1][2]) / A_det;
        A_inv[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) / A_det;
        A_inv[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) / A_det;
        A_inv[1][0] = -(A[1][0] * A[2][2] - A[1][2] * A[2][0]) / A_det;
        A_inv[1][1] = (A[0][0] * A[2][2] - A[0][2] * A[2][0]) / A_det;
        A_inv[1][2] = -(A[0][0] * A[1][2] - A[0][2] * A[1][0]) / A_det;
        A_inv[2][0] = (A[1][0] * A[2][1] - A[1][1] * A[2][0]) / A_det;
        A_inv[2][1] = -(A[0][0] * A[2][1] - A[0][1] * A[2][0]) / A_det;
        A_inv[2][2] = (A[0][0] * A[1][1] - A[0][1] * A[1][0]) / A_det;

        double v1[N], v2[N], v3[N];
        matrix_multiply(A_inv, B, v1);
        matrix_multiply(A_inv, C, v2);
        v3[0] = 1;
        v3[1] = v3[2] = 0;

        matrix_transpose(v, v);
        for (i = 0; i < N; i++) {
            v[0][i] = v1[i];
            v[1][i] = v2[i];
            v[2][i] = v3[i];
        }
    } else if (B_det != 0) {
        double B_inv[N][N];
        B_inv[0][0] = (B[1][1] * B[2][2] - B[2][1] * B[1][2]) / B_det;
        B_inv[0][1] = -(B[0][1] * B[2][2] - B[0][2] * B[2][1]) / B_det;
        B_inv[0][2] = (B[0][1] * B[1][2] - B[0][2] * B[1][1]) / B_det;
        B_inv[1][0] = -(B[1][0] * B[2][2] - B[1][2] * B[2][0]) / B_det;
        B_inv[1][1] = (B[0][0] * B[2][2] - B[0][2] * B[2][0]) / B_det;
        B_inv[1][2] = -(B[0][0] * B[1][2] - B[0][2] * B[1][0]) / B_det;
        B_inv[2][0] = (B[1][0] * B[2][1] - B[1][1] * B[2][0]) / B_det;
        B_inv[2][1] = -(B[0][0] * B[2][1] - B[0][1] * B[2][0]) / B_det;
        B_inv[2][2] = (B[0][0] * B[1][1] - B[0][1] * B[1][0]) / B_det;

        double v1[N], v2[N], v3[N];
        matrix_multiply(B_inv, C, v1);
        v2[0] = 1;
        v2[1] = v2[2] = 0;
        v3[0] = 0;
        v3[1] = 1;
        v3[2] = 0;

        matrix_transpose(v, v);
        for (i = 0; i < N; i++) {
            v[0][i] = v1[i];
            v[1][i] = v2[i];
            v[2][i] = v3[i];
        }
    } else if (C_det != 0) {
        double C_inv[N][N];
        C_inv[0][0] = (C[1][1] * C[2][2] - C[2][1] * C[1][2]) / C_det;
        C_inv[0][1] = -(C[0][1] * C[2][2] - C[0][2] * C[2][1]) / C_det;
        C_inv[0][2] = (C[0][1] * C[1][2] - C[0][2] * C[1][1]) / C_det;
        C_inv[1][0] = -(C[1][0] * C[2][2] - C[1][2] * C[2][0]) / C_det;
        C_inv[1][1] = (C[0][0] * C[2][2] - C[0][2] * C[2][0]) / C_det;
        C_inv[1][2] = -(C[0][0] * C[1][2] - C[0][2] * C[1][0]) / C_det;
        C_inv[2][0] = (C[1][0] * C[2][1] - C[1][1] * C[2][0]) / C_det;
        C_inv[2][1] = -(C[0][0] * C[2][1] - C[0][1] * C[2][0]) / C_det;
        C_inv[2][2] = (C[0][0] * C[1][1] - C[0][1] * C[1][0]) / C_det;

        double v1[N], v2[N], v3[N];
        v1[0] = 0;
        v1[1] = 0;
        v1[2] = 1;
        v2[0] = 0;
        v2[1] = 1;
        v2[2] = 0;
        v3[0] = 1;
        v3[1] = 0;
        v3[2] = 0;

        matrix_transpose(v, v);
        for (i = 0; i < N; i++) {
            v[0][i] = v1[i];
            v[1][i] = v2[i];
            v[2][i] = v3[i];
        }
    } else {
        printf("Error: the matrix is singular!\n");
    }
}

int main() {
    double a[N][N] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    double lambda[N];
    double v[N][N];

    eigenvalues(a, lambda);
    eigenvectors(a, lambda, v);

    printf("Eigenvalues:\n");
    printf("%f %f %f\n", lambda[0], lambda[1], lambda[2]);

    printf("Eigenvectors:\n");
    matrix_print(v);

    return 0;
}

这个示例程序中，首先定义了一个 `matrix_multiply()` 函数来实现矩阵乘法，一个 `matrix_transpose()` 函数来实现矩阵转置，一个 `matrix_print()` 函数来打印矩阵。然后定义了一个 `eigenvalues()` 函数来求解特征值，一个 `eigenvectors()` 函数来求解特征向量。最后在 `main()` 函数中给出了一个三阶矩阵作为例子，并打印出特征值和特征向量。

值得注意的是，这个示例程序中的特征值和特征向量的计算使用了三种不同的方法，分别适用于三种不同的情况。具体来说，当矩阵 A 的行列式不为 0 时，我们可以通过求解矩阵 (A - λI) 的逆矩阵来求解特征向量；当矩阵 B 的行列式不为 0 时，我们可以通过求解矩阵 (B - λI) 的逆矩阵来求解特征向量；当矩阵 C 的行列式不为 0 时，我们可以通过求解矩阵 (C - λI) 的逆矩阵来求解特征向量；当矩阵 A、B 和 C 的行列式均为 0 时，矩阵是奇异矩阵，无法求解特征向量。