高斯消除:逆 4x4 矩阵

Gaussian Elimination: Inverse 4x4 Matrix

本文关键字:4x4 矩阵 高斯消      更新时间:2023-10-16

问题

下面会给我一个变量"cur"的单位矩阵,

并试图给我变量"tmp"的逆矩阵,但会失败。谁能告诉我我做错了什么?

附言矩阵是列主矩阵。

矩阵输入

1 | 0

| 0 | 0

0 | -

4.37114e-08 | 1 | 0

0 | -1 | -

4.37114e-08 | 0

0 | 0 | 0 | 1

矩阵输出

1 | 0

| 0 | 0

0 | 0 | 0

| 0 | 0

0 | 1 | 0 | 0

0 | 0 | 0 | 1

期望的输出

1 | 0

| 0 | 0

0 | -

4.3711399999999916e-8 | -0.9999999999999981 | 0

0 | 0.999999999999998 | -

4.3711399999999916e-8 | 0

0 | 0 | 0 | 1

#include <iostream>
template <typename T>
class Matrix
{
public:
    T matrix[4][4];
    Matrix(T matrix[4][4])
        : matrix()
    {
        for (unsigned int y = 0; y < 4; ++y)
        {
            for (unsigned int x = 0; x < 4; ++x)
            {
                this->matrix[y][x] = matrix[y][x];
            }
        }
    }
    Matrix()
        : matrix()
    {
        T zero = static_cast<T>(0);
        T one = static_cast<T>(1);
        matrix[0][0] = one;
        matrix[1][0] = zero;
        matrix[2][0] = zero;
        matrix[3][0] = zero;
        matrix[0][1] = zero;
        matrix[1][1] = one;
        matrix[2][1] = zero;
        matrix[3][1] = zero;
        matrix[0][2] = zero;
        matrix[1][2] = zero;
        matrix[2][2] = one;
        matrix[3][2] = zero;
        matrix[0][3] = zero;
        matrix[1][3] = zero;
        matrix[2][3] = zero;
        matrix[3][3] = one;
    }
    Matrix<T> GetInverse() const
    {
        T zero = static_cast<T>(0);
        T one = static_cast<T>(1);
        Matrix<T> tmp;
        Matrix<T> cur;
        for (unsigned int y = 0; y < 4; ++y)
        {
            for (unsigned int x = 0; x < 4; ++x)
            {
                cur.matrix[y][x] = matrix[y][x];
            }
        }
        cur.Print();
        for (unsigned int x = 0; x < 4; x++)
        {
            if (cur.matrix[x][x] != zero)
            {
                T denominator = cur.matrix[x][x];
                for (unsigned int a = x; a < 4; ++a)
                {
                    cur.matrix[x][a] = cur.matrix[x][a] / denominator;
                }
                for (unsigned int a = 0; a < 4; ++a)
                {
                    tmp.matrix[x][a] = tmp.matrix[x][a] / denominator;
                }
            }
            for (unsigned int y = 0; y < 4; ++y)
            {
                if (y != x && cur.matrix[y][x] != 0)
                {
                    T difference = cur.matrix[y][x];
                    for (unsigned int a = x; a < 4; ++a)
                    {
                        cur.matrix[y][a] = (cur.matrix[y][a] - difference) * cur.matrix[x][a];
                    }
                    for (unsigned int a = 0; a < 4; ++a)
                    {
                        tmp.matrix[y][a] = (tmp.matrix[y][a] - difference) * tmp.matrix[x][a];
                    }
                }
            }
        }
        cur.Print();
        tmp.Print();
        return tmp;
    }
    void Print()
    {
        for (unsigned int y = 0; y < 4; ++y)
        {
            for (unsigned int x = 0; x < 4; ++x)
            {
                std::cout << matrix[y][x];
                if (x < 3)
                    std::cout << " | ";
            }
            std::cout << std::endl;
        }
        std::cout << std::endl;
    }
}
int main()
{
    float matrix[4][4];
    matrix[0][0] = 1.0f;
    matrix[0][1] = 0.0f;
    matrix[0][2] = 0.0f;
    matrix[0][3] = 0.0f;
    matrix[1][0] = 0.0f;
    matrix[1][1] = -4.37114e-08;
    matrix[1][2] = 1.0f;
    matrix[1][3] = 0.0f;
    matrix[2][0] = 0.0f;
    matrix[2][1] = -1.0f;
    matrix[2][2] = -4.37114e-08;
    matrix[2][3] = 0.0f;
    matrix[3][0] = 0.0f;
    matrix[3][1] = 0.0f;
    matrix[3][2] = 0.0f;
    matrix[3][3] = 1.0f;
    Matrix<float> inverseMatrix(matrix);
    inverseMatrix.GetInverse();
    return 0;
}

首先,你的类不是一个合适的类。在尝试找出算法或操作错误之前,您可以将其转换为更合理的代码,以便于调试。

现在,您的Matrix类并没有真正定义对象,而只是一组要在 2D 数组上使用的操作。一个好的对象具有属性和数据,然后具有使用该数据的特殊函数。我会在 https://en.wikipedia.org/wiki/Object-oriented_programming 阅读一些介绍

你想要的是类似的东西

// Matrix.h
template <typename T>
class Matrix {
public:
    // exactly how you choose to construct the object is up to you
    // i.e. if you want to pass a 2d array as a parameter or something
    Matrix(T inputMatrix[4][4]) {
        for (int column = 0; column < 4; ++column) {
            for (int row = 0; row < 4; ++row) {
                matrix[column][row] = inputMatrix[column][row];
            }
        }
    } 
    Matrix reduce() {
        // will return a new Matrix object
        . . .
    }
    Matrix inverse() {
        // probably will end up calling reduce()
        . . .
    }
    void print() {
        . . .
    }
    // any other operations you want to do on matrices
private:
    T matrix[4][4];
};

因此,当您想进行操作时,它非常简单明了:

#include "Matrix.h"
int main() {
    float matrix[4][4] = {
        { 1.0f, 0.0f, 0.0f, 0.0f },
        { 0.0f, -4.37114e-08, 1.0f, 0.0f },
        { 0.0f, -1.0f, -4.37114e-08, 0.0f },
        { 0.0f, 0.0f, 0.0f, 0.0f };
    Matrix<float> A = Matrix<float>(matrix);
    Matrix<float> B = A.inverse();
    B.print();
    return 0;
}