一般行列式求值

最新推荐文章于 2024-09-19 15:37:02 发布

crushed-sw

最新推荐文章于 2024-09-19 15:37:02 发布

阅读量412

点赞数 1

文章标签：线性代数矩阵算法

本文链接：https://blog.csdn.net/qq_50952332/article/details/125020802

版权

求一般行列式的值

求行列式的值

让我们来复习一下行列式的值是怎么求的，有对角线法，代数余子式法，等价转换法，逆序数法等方法。但有的计算麻烦，有的只适合特殊的行列式使用，没有普适性，因此这里只介绍一种具有普适性的方法以供参考。该法便是全选主元高斯消去法。

高斯消去法就是对方阵 $A$ 进行一系列的变化，使之成为上三角矩阵或下三角矩阵，其对角线的各个元素乘积即为行列式的值。

变换过程如下：

对于 $\dots, n-2$ 作以下变换：
$a_{ij}- \frac{a_{ik}a_{kj}}{a_{kk}} \Rightarrow a_{ij} \qquad i, j = k+1, \dots, n-1\\$
为了保证计算时的稳定性，在实际变换过程中使用全选主元的方法，全选主元的意思就是将绝对值最大的元素交换到主对角线上，全选主元的过程如下所示：

对于矩阵求逆过程中的第 $k$ 步，首先，在 $a_{kk}$ 的右下方（包括 $a_{kk}$ 自己）的 $n - k + 1$ 阶子矩阵中选取绝对值最大的元素，并将该元素所在的行号记录在 $I S (k)$ 中，而列号记录在 $J S (k)$ 中，然后通过行交换与列交换将该绝对值最大的元素交换到主对角线 $a_{kk}$ 的位置上，即做以下两步：
$\Leftrightarrow A(IS(k),l),l=1,2,...,n\\ A(l,k) \Leftrightarrow A(l,JS(k)),l=1,2,...,n\\$
经过全选主元过后，矩阵求逆的过程是稳定的。但最后需要对结果进行恢复。恢复过程如下：
$\Leftrightarrow A(JS(k),l),l=1,2,...,n\\ A(l,k) \Leftrightarrow A(l,IS(k)),l=1,2,...,n\\$

我们对其接受进行设计：double det(double* a, int n)

接口与参数说明：

形参与函数类型	参数意义
`double* a`	存放矩阵 $A$ 的元素，返回时会被破坏
`const int n`	存放矩阵的阶数
`double det()`	返回行列式的值

这时便可以轻松写出该程序：

#include <math.h>

template <int N>
double det(double a[N][N]) {
    int i, j, k;
    int is, js;
    double f, detValue, q, d;

    f = 1.0;
    detValue = 1.0;

    for(k = 0; k != N-1; ++k) {
        q = 0.0;
        // 遍历右下子矩阵，找出最大元素
        for(i = k; i != N; ++i) {
            for(j = k; j != N; ++j) {
                d = fabs(a[i][j]);
                if(d > q) {
                    q = d;
                    is = i;
                    js = j;
                }
            }
        }

        // 全为0的矩阵
        if(q + 1.0 == 1.0) {
            detValue = 0.0;
            return detValue;
        }

        // 最大元素不在主对角线，
        // 进行逐行变换
        if(is != k) {
            f = -f;
            for(j = k; j != N; ++j) {
                d = a[k][j];
                a[k][j] = a[is][j];
                a[is][j] = d;
            }
        }

        // 逐列变换
        if(js != k) {
            f = -f;
            for(i = k; i != N; ++i) {
                d = a[k][i];
                a[k][i] = a[js][i];
                a[js][i] = d;
            }
        }

        // 将变换后的主对角线元素乘到结果中
        detValue *= a[k][k];

        // 行变换消元
        for(i = k+1; i != N; ++i) {
            d = a[i][k] / a[k][k];
            for(j = k+1; j != N; ++j) {
                a[i][j] -= d * a[k][j];
            }
        }
    }

    // 将最后一行的元素乘到结果中
    detValue = f * detValue * a[N-1][N-1];
    return detValue;
}

在这里来测试一下这个程序：

$\begin{bmatrix} 1 & 2 & 3 & 4\\ 5 & 6 & 7 & 8\\ 9 & 10 & 11 & 12\\ 13 & 14 & 15 & 16\\ \end{bmatrix}$

$\begin{bmatrix} 3 & -3 & -2 & 4\\ 5 & -5 & 1 & 8\\ 11 & 8 & 5 & -7\\ 5 & -1 & -3 & -1\\ \end{bmatrix}$

运行程序如下：

#include "det.hpp"

int main() {
    double a[4][4] = {{1, 2, 3, 4},
                      {5, 6, 7, 8},
                      {9, 10, 11, 12},
                      {13, 14, 15, 16}};

    double b[4][4] = {{3, -3, -2, 4},
                      {5, -5, 1, 8},
                      {11, 8, 5, -7},
                      {5, -1, -3, -1}};

    std::cout << "det(a)=" << det<4>(a) << std::endl;
    std::cout << "det(b)=" << det<4>(b) << std::endl;

    return 0;
}

运行结果如下所示：