Octave高斯消元法函数报数组越界错误res(2): out of bound 1求助

问题现象 #

在Octave中运行高斯消元法代码时，gauss函数末尾打印结果触发错误：

error: res(2): out of bound 1 (dimensions are 1x1)

原始代码 #

main.m #

n9=9;
n10=10;
A1 = [61, -82, 35, -48, -40, -43, -50, -1, -92;
      -39, -99, -75, -75, 26, -1, 42, 21, 10;
      85, -88, 76, 48, -8, 90, 94, -86, -24;
      28, 8, -31, -29, 94, 98, -95, 96, 50;
      86, -37, 12, 14, -99, -26, 33, 17, 88;
      24, -30, -51, -57, 14, -53, 58, 21, -92;
      7, -9, -43, 69, -82, 56, -30, 100, -44;
      24, 37, 98, -77, -55, 16, 41, 12, 46;
      81, 79, 8, -70, 65, -76, -39, -82, -54]

A2 = [268, 74, -71, 11, -104, 16, -33, 64, 7, 11;
      74, 375, -118, 230, -51, -212, 180, 204, 163, -35;
      -71, -118, 400, -122, -81, 72, -201, 110, 20, 333;
      11, 230, -122, 239, -47, -40, 189, 188, 57, -69;
      -104, -51, -81, -47, 325, -53, -13, -222, -27, -118;
      16, -212, 72, -40, -53, 335, -64, -60, -110, 73;
      -33, 180, -201, 189, -13, -64, 325, 47, -8, -84;
      64, 204, 110, 188, -222, -60, 47, 446, 123, 123;
      7, 163, 20, 57, -27, -110, -8, 123, 379, 65;
      11, -35, 333, -69, -118, 73, -84, 123, 65, 378]

A3 = [-327, 106, 0, 0, 0, 0, 0, 0, 0, 0;
      -104, -512, 92, 0, 0, 0, 0, 0, 0, 0;
      0, -5, 56, 6, 0, 0, 0, 0, 0, 0;
      0, 0, 42, 501, 26, 0, 0, 0, 0, 0;
      0, 0, 0, 2, 167, -70, 0, 0, 0, 0;
      0, 0, 0, 0, 218, 439, -64, 0, 0, 0;
      0, 0, 0, 0, 0, -145, 943, 365, 0, 0;
      0, 0, 0, 0, 0, 0, -8, -91, -26, 0;
      0, 0, 0, 0, 0, 0, 0, -143, -749, 169;
      0, 0, 0, 0, 0, 0, 0, 0, 219, 566]
b1 = [67, 94, -83, 76, -93, -79, -18, 94, -77]
b2 = [420, 19, 384, 83, 108, 264, 80, 14, -215, 394]
b3 = [-981, 687, 845, 542, -915, -244, 409, 459, -552, -462]

gauss(A1, b1, n9)

gauss.m #

function result = gauss (a, b, n)
  res = [9];
  temp = 0;
  for k = 1:n,
    for i = k+1:n,
      for j = k:n,
        if j == k,
          temp = a(i,j) ./ a(k,j);
          a(i,j) = a(i,j) - a(k,j) .* a(i,j) ./ a(k,j);
          b(i) = b(i) - b(k) * temp;
        else
          a(i,j) = a(i,j) - temp .* a(k,j);
        endif;
      end;
    end;
  end;
  for i = n:1,
    for j = n-1:0,
      if a(i, j) != 0,
        if i == j,
           res = b(j) ./ a(i,j);
        else
           b(i) = b(i) - a(i,j) .* res(j);
        endif;
      endif;
    end;
  end;

  for i = 1:9,
    printf("%f", res(i))
  end;

endfunction

问题分析 #

1. res数组初始化错误：res = [9];创建的是1x1的数组，无法存储n个解，应初始化为长度为n的数组。
2. 回代循环索引错误：
   • for i = n:1步长默认是1，n到1无法生成有效循环序列，需改为for i = n:-1:1（从n倒序到1）。
   • j = n-1:0会生成0索引，而Octave数组索引从1开始，且循环方向错误，应改为for j = i+1:n（遍历i之后的列）。
3. res赋值错误：res = b(j) ./ a(i,j);会将res覆盖为单个值，而非赋值到对应索引位置，应改为res(i) = b(i) ./ a(i,i);。
4. 打印循环硬编码：for i = 1:9固定了循环次数，应使用参数n适配不同维度的矩阵。
5. 消元步骤冗余：a(i,j) = a(i,j) - a(k,j) .* a(i,j) ./ a(k,j);等价于a(i,j) = 0;，可直接简化。

修正后的gauss.m #

function result = gauss (a, b, n)
  res = zeros(1, n); % 初始化长度为n的结果数组
  temp = 0;
  
  % 前向消元
  for k = 1:n,
    for i = k+1:n,
      temp = a(i,k) / a(k,k); % 计算消元因子
      a(i,k:n) = a(i,k:n) - temp * a(k,k:n); % 整行消元
      b(i) = b(i) - temp * b(k);
    end;
  end;
  
  % 回代求解
  for i = n:-1:1,
    res(i) = b(i);
    for j = i+1:n,
      res(i) = res(i) - a(i,j) * res(j);
    end;
    res(i) = res(i) / a(i,i);
  end;
  
  % 打印结果
  for i = 1:n,
    printf("x%d = %.6f\n", i, res(i));
  end;
  
  result = res; % 返回结果数组
endfunction

修正说明 #

• 初始化res为zeros(1, n)，确保数组长度与解的数量匹配。
• 简化前向消元步骤，直接对整行进行操作，避免冗余计算。
• 修正回代循环的索引方向，从最后一行开始倒序计算，正确累加已求出的解。
• 打印循环使用参数n，适配9阶、10阶等不同维度的矩阵。
• 最后将结果赋值给result，确保函数能返回求解结果。

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