使用部分旋转实现高斯消除

Question

我正在编写一个程序来在 MATLAB 中实现具有部分旋转的高斯消去法。我创建了一个整数数组来存储行的交换，而不是直接交换行。

但是，我无法得到正确的结果，也无法找出问题所在。我应该如何修改我的代码以获得正确答案？

A=[1, -2, 1; 2, 1, -3; 4, -7, 1];
b=[0, 5, 1];
b=b';

function x = GE(A,b)

[m,n]= size(A); 
if m ~= n
  disp('Not a square');
end

for p=1:n
   array(p)=p;
endfor

A = [A,b];      

%elimination
for i = 1:n-1
    pivot = i;
    %select pivot
    for j = i+1:n
        if  abs(A(array(i),i)) < abs(A(array(i),i)) %row interchange
            temp = array(i);
            array(i) = array(j);
            array(j) = temp;
        end
    end

    while (pivot <= n && A(pivot,i)== 0)
        pivot = pivot+1;
    end
    if pivot > n
        disp('No unique solution');
        break
    else
        if pivot > i
            tem = array(i);
            array(i) = pivot
            pivot= tem;
        end
    end

    for j = i+1:n
        m = -A(array(j),i)/A(array(i),i);
        for k = i+1:n+1 
            A(array(j),k) = A(array(j),k) + m*A(array(i),k);
        end
    end
end

if A(n,n) == 0
    disp('No unique solution');
    return
end

%backward substitution
x(n) = A(array(n),n+1)/A(array(n),n);
for i = n - 1:-1:1
    sum = 0;
    for j = i+1:n
        sum = sum + A(array(i),j)*x(j);
    end
    x(i) = (A(array(i),n+1) - sum)/A(array(i),i);
end

endfunction

Answer 1

你的错误很简单。你没有正确旋转 - 特别是在你的 if 声明中：

for j = i+1:n
    if  abs(A(array(i),i)) < abs(A(array(i),i)) %row interchange <-------
        temp = array(i);
        array(i) = array(j);
        array(j) = temp;

    end
end

您检查以哪个系数为基准的检查不正确，因为您在需要使用 j 时使用索引 i 沿列检查。因此，只需更改布尔条件的另一端，使其使用 j:

for j = i+1:n
    if  abs(A(array(i),i)) < abs(A(array(j),j)) %// CHANGE HERE
        temp = array(i);
        array(i) = array(j);
        array(j) = temp;

    end
end

以你的例子A和b，我得到的答案是x = [2 1 0]，这与MATLAB的ldivide操作一致：x = A \ b.

使用部分旋转实现高斯消除

Implementing Gaussian elimination with partial pivoting

matlab

gaussian