(C++)n阶方阵求逆

一、实验目的、内容

输入是一个$n$（$n<256$）阶方阵$A$，输出是它的逆矩阵，再将得到的逆矩阵与原来的矩阵相乘，验证其结果是单位矩阵。

二、实验程序设计及结构

1.需求分析

变量

二级指针a(double**)存储原矩阵的首地址，b(double**)存储逆矩阵的首地址，c(double**)存储乘积的首地址；数n(unsigned char)存储阶数；循环变量i(unsigned char)，j(unsigned char)。

函数

求逆函数double** inv(const double** a, const unsigned char n)，求积函数double** pro(const double** a, const double** b, const unsigned char n)，主函数int main()。

2.设计结构或流程图

1. 输入原矩阵及矩阵的阶数。
2. 求逆矩阵并输出。（调用求逆函数）
   1. 在堆区开辟两个二维数组，其中一个初始化为原矩阵，另一个初始化为单位矩阵。
   2. 从首行开始逐行进行初等行变换，将原矩阵化为阶梯型，对单位矩阵进行同等操作。
   3. 从最后一行开始逐行消元素，直到将原矩阵化为单位矩阵为止。
3. 分别对原矩阵和逆矩阵左乘和右乘进行验证并输出。（调用求积函数）
4. 释放堆区数据并退出。

三、设计过程

#include <iostream>
#define abs(a) (a > 0 ? a : -a)
using namespace std;
// 求逆函数
double **inv(double **a, unsigned char n)
{
    unsigned char i, j, k = -1;
    double **t = new double *[n], **b = new double *[n]; // t作为逆矩阵，b作为原矩阵的副本
    // 拷贝a并把t初始化为单位矩阵
    for (i = 0; i < n; ++i)
    {
        t[i] = new double[n];
        b[i] = new double[n];
        for (j = 0; j < n; ++j)
        {
            b[i][j] = a[i][j];
            if (i == j)
                t[i][j] = 1.;
            else
                t[i][j] = 0.;
        }
    }
    double x;
    double *w;
    while (++k < n) // 化原矩阵为阶梯型
    {
        // 第k行第k列为零，在第k列中寻找不为零的元素
        if (abs(b[k][k]) < 1e-5)
        {
            i = k;
            // 由于浮点数为二进制指数表达，在计算过程中可能出现偏差
            // 所以将1e-5到-1e-5之间的数据近似地视为0
            do
                if (++i == n)
                    return nullptr;
            while (abs(b[i][k]) < 1e-5); // 如果遍历完了整行，证明行列式为0，逆矩阵不存在
            w = b[i];
            b[i] = b[k];
            b[k] = w; // 交换找到的行
            w = t[i];
            t[i] = t[k];
            t[k] = w; // 对逆矩阵做同样的初等行变换
        }
        for (i = k + 1; i < n; ++i)
        {
            x = b[i][k] / b[k][k]; // 记录比例系数
            for (j = k + 1; j < n; ++j)
                b[i][j] -= x * b[k][j]; // 对原矩阵做初等行变换
            for (j = 0; j <= k; ++j)
                t[i][j] -= x * t[k][j]; // 对逆矩阵做同样的初等行变换
        }
        x = b[k][k];
        // 将b[k][k]变为1，由于b[k][k]在接下来的程序里无用，故i从k+1开始
        for (i = k + 1; i < n; ++i)
            b[k][i] /= x;
        for (i = 0; i <= k; ++i)
            t[k][i] /= x; // 对逆矩阵做同样的初等行变换
    }
    do // 化为单位矩阵
    {
        i = --k;
        do
        {
            x = b[--i][k];
            for (j = 0; j < n; ++j)
                t[i][j] -= x * t[k][j];
        } while (i > 0);
        delete[] b[k]; // 防止内存泄漏
    } while (k > 1);
    delete[] *b;
    delete[] b;
    return t;
}
double **pro(double **a, double **b, unsigned char n) // 方阵乘法
{
    double **p = new double *[n], *q;
    for (unsigned char i = 0; i < n; ++i)
    {
        q = p[i] = new double[n];
        for (unsigned char j = 0; j < n; ++j)
        {
            *q = *a[i] * b[0][j];
            for (unsigned char k = 1; k < n; ++k)
                *q += a[i][k] * b[k][j];
            ++q;
        }
    }
    return p;
}
int main()
{
    unsigned short n; // 为了输入方便，故采用unsigned short
    cout << "请输入方阵的阶数:\n";
    cin >> n;
    cout << "请输入方阵:\n";
    double **a = new double *[n];
    // 由于a是指针类型而非数组类型，故不能用范围for循环
    for (unsigned char i = 0; i < n; ++i)
    {
        a[i] = new double[n];
        for (unsigned char j = 0; j < n; ++j)
            cin >> a[i][j];
    }
    double **b = inv(a, n);
    if (!b)
    {
        cout << "逆矩阵不存在!\n";
        system("pause");
        return 0;
    }
    cout << "逆矩阵为:\n";
    for (unsigned char i = 0; i < n; ++i)
    {
        for (unsigned char j = 0; j < n; ++j)
            cout << b[i][j] << '\t';
        cout << endl;
    }
    cout << "原矩阵右乘逆矩阵为:\n";
    double **c = pro(a, b, n);
    for (unsigned char i = 0; i < n; ++i)
    {
        for (unsigned char j = 0; j < n; ++j)
            cout << c[i][j] << '\t';
        delete[] c[i]; // 防止内存泄漏
        cout << endl;
    }
    delete[] c;
    cout << "原矩阵左乘逆矩阵为:\n";
    c = pro(b, a, n);
    for (unsigned char i = 0; i < n; ++i)
    {
        for (unsigned char j = 0; j < n; ++j)
            cout << c[i][j] << '\t';
        delete[] c[i]; // 防止内存泄漏
        delete[] a[i];
        delete[] b[i];
        cout << endl;
    }
    delete[] c;
    delete[] a;
    delete[] b;
    system("pause");
    return 0;
}

四、测试分析

第一组

第二组

实验中出现的bug及解决方案

bug	解决方案
未考虑行列式为0的情况	返回空指针
由于浮点数计算的不精确，容易出现大数	将-1e-5至1e-5之间的数据近似视为0

五、设计的特点和结果

采用初等变换法求逆矩阵，较伴随矩阵法较复杂，但避免了递归，效率更高。

结果：求得逆矩阵。但由于浮点数运算的偏差，容易出现极小数，往后可以通过编写一个有理数类加以解决。