Octave中Gaussian elimination执行后出现NaN问题求助

代码里的问题与修正 #

1. 循环变量笔误 #

在Forward Elimination阶段，i=k+i:n是明显的书写错误，必须改成i=k+1:n，否则会导致循环变量i迅速超出矩阵范围，引发计算错误。

2. 消元因子覆盖原始数据 #

原代码直接修改a(i,k)为消元因子，后续行的消元计算会误用这个修改后的值，逻辑完全错误。正确做法是先把消元因子存到临时变量中，再用该变量更新矩阵和b向量，避免覆盖主元列的原始数据。

3. 矩阵本身是奇异矩阵 #

你使用的矩阵a是奇异矩阵（任意相邻两行的差值都是[1 1 1 1 1]，行向量线性相关，行列式为0），这类矩阵不存在唯一解。即使代码修正，也会出现主元为0的情况，导致除以0得到NaN。若必须使用该矩阵，需明确它要么无解、要么有无穷多解；若仅为测试代码，建议更换非奇异矩阵。

修正后的完整代码 #

clc; clear all; close all;

a = [1 2 3 4 5;
     2 3 4 5 6;
     3 4 5 6 7;
     4 5 6 7 8;
     5 6 7 8 9;]

b = [1; 2; 3; 4; 5;]
n = size(a,1);

% Gaussian Elimination
for k=1:n-1
  % 提前检查主元是否接近0，防止除以0
  if abs(a(k,k)) < 1e-10
    error('矩阵奇异，无法进行高斯消元');
  endif
  for i=k+1:n
    % 用临时变量存储消元因子，避免覆盖原始数据
    factor = a(i,k)/a(k,k);
    a(i,k) = factor;
    for j=k+1:n
      a(i,j)= a(i,j) - factor*a(k,j);
    endfor
  endfor
endfor
disp("a-----------------------")
disp(a);

% Forward Elimination
for k=1:n-1
   for i=k+1:n  % 修正循环变量
     b(i) = b(i) - a(i,k)*b(k);
   endfor
endfor

% Backward solve
x = zeros(n,1);
for i=n:-1:1
   s=b(i);
   for j=i+1:n
     s=s-a(i,j)*x(j);
   endfor
   % 回代时检查主元，避免除以0
   if abs(a(i,i)) < 1e-10
     error('主元为0，无法回代求解');
   endif
   x(i)= s/a(i,i);
endfor
disp(x);

额外测试提示 #

运行修正后的代码会弹出“矩阵奇异”的错误，这是原矩阵的固有问题。若要验证高斯消元逻辑的正确性，可以更换一个非奇异矩阵，比如：

a = [2 1 -1; -3 -1 2; -2 1 2];
b = [8; -11; -3];

这个矩阵的唯一解为x = [2; 3; -1]，能正常运行得到结果。

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