基于Gauss-Jordan法的Mathematica矩阵求逆程序错误排查

你实现的Gauss-Jordan法求逆矩阵代码存在多处逻辑错误，导致结果不是原矩阵的逆，但因行列式的乘积性质（可逆矩阵的逆矩阵行列式是原矩阵行列式的倒数，乘积自然为1），所以出现了行列式乘积为1但矩阵不是逆的矛盾情况。具体错误及修正方案如下：

核心错误点 #

1. 单位矩阵初始化逻辑完全错误 #

你初始化存储逆矩阵的mat2时，循环内的赋值逻辑混乱：

For[i = 1, i <= m, i++,
  For[j = 1, j <= n, 
   j++, {Subscript[c, i, i] = 1, Subscript[c, i, j] = 0, 
    Subscript[c, m, n] = 1}]];

这段代码在每次j循环时都会把Subscript[c,i,j]设为0（包括j=i的对角线位置），最后还强行给Subscript[c,m,n]赋值1，导致mat2根本不是单位矩阵，后续所有变换都失去了正确的基准。

2. 上三角消元的索引逻辑颠倒 #

最后一步消去主元上方元素的循环中，错误地操作了行索引而非目标行：

For[t = k, t <= n, t++,
 Subscript[b, t, l] = Subscript[b, t, l] - coe3*Subscript[b, t, k];
 Subscript[c, t, l] = 
  Subscript[c, t, l] - coe3*Subscript[c, t, k]]

Gauss-Jordan的最后一步应该是用已归一化的第l行，消去第k行的第l列元素，即对第k行进行整体操作，而非遍历t行修改第l列。

3. 语法错误与冗余操作 #

• Mathematica中没有小写的continue，正确语句是Continue[]；
• 手动给下标变量赋值随机数的方式冗余，直接生成矩阵更简洁；
• 未限制m和n必须相等（求逆仅适用于方阵）。

修正后的代码 #

n = Input["请输入方阵阶数n："];
(* 初始化原矩阵和单位矩阵 *)
mat = RandomInteger[{1, 3}, {n, n}];
mat2 = IdentityMatrix[n];
Print["原矩阵："];
Print[mat // MatrixForm];

(* 将矩阵元素映射到下标变量，方便分步操作 *)
Do[Subscript[b, i, j] = mat[[i, j]], {i, n}, {j, n}];
Do[Subscript[c, i, j] = mat2[[i, j]], {i, n}, {j, n}];

(* Gauss-Jordan一步消元：同时处理主元上下的元素 *)
For[k = 1, k <= n, k++,
  (* 主元选择：避免除以0，增强鲁棒性 *)
  maxRow = First[Ordering[Abs[Table[Subscript[b, i, k], {i, k, n}]], -1]] + k - 1;
  If[maxRow != k,
    (* 交换原矩阵的行 *)
    Do[temp = Subscript[b, k, t]; Subscript[b, k, t] = Subscript[b, maxRow, t]; Subscript[b, maxRow, t] = temp, {t, n}];
    (* 交换逆矩阵的行 *)
    Do[temp = Subscript[c, k, t]; Subscript[c, k, t] = Subscript[c, maxRow, t]; Subscript[c, maxRow, t] = temp, {t, n}];
  ];
  (* 归一化主元所在行 *)
  coe = 1/Subscript[b, k, k];
  Do[Subscript[b, k, t] = Subscript[b, k, t] * coe; Subscript[c, k, t] = Subscript[c, k, t] * coe, {t, n}];
  (* 消去当前主元列的其他所有元素 *)
  Do[If[i != k,
    coe2 = Subscript[b, i, k];
    Do[Subscript[b, i, t] = Subscript[b, i, t] - coe2 * Subscript[b, k, t];
       Subscript[c, i, t] = Subscript[c, i, t] - coe2 * Subscript[c, k, t], {t, n}];
  ], {i, n}];
];

(* 将下标变量转换回矩阵 *)
invMat = Table[Subscript[c, i, j], {i, n}, {j, n}];
Print["原矩阵的逆："];
Print[invMat // MatrixForm];
Print["验证：原矩阵×逆矩阵（应为单位矩阵）："];
Print[mat . invMat // MatrixForm];
Print["行列式乘积：", Det[mat] * Det[invMat]];
Quit[];

修正说明 #

1. 用IdentityMatrix直接生成单位矩阵，避免手动循环的错误；
2. 加入主元选择逻辑，解决原矩阵主元为0的异常情况；
3. 采用Gauss-Jordan一步消元法，同时消去主元上下的元素，逻辑更简洁；
4. 修正了循环索引的错误，确保对目标行进行正确操作；
5. 增加验证步骤，直接输出原矩阵与逆矩阵的乘积，确认结果是否为单位矩阵。

内容的提问来源于stack exchange，提问作者yasirbagci