Gauss-Seidel 方法不适用于大型稀疏数组？

Question

我再次遇到了 Matlab 中的 Gauss-Seidel 方法的问题。这里是：

function [x] = ex1_3(A,b)

format long 

sizeA=size(A,1);

x=zeros(sizeA,1);  

%Just a check for the conditions of the Gauss-Seidel Method (if it has dominant diagonal)

for i=1:sizeA
    sum=0;
    for j=1:sizeA
        if i~=j 
            sum=sum+abs(A(i,j));
        end
    end
    if abs(A(i,i))<sum
        fprintf('\nGauss-Seidel''s conditions not met!\n');
        return     
    end
end

%Actual Gauss-Seidel Method

max_temp=10^(-6); %Pass first iteration
while max_temp>(0.5*10^(-6))
    xprevious=x;
    for i=1:sizeA
        x(i,1)=b(i,1);
        for j=1:sizeA
            if i~=j
              x(i,1)=x(i,1)-A(i,j)*x(j,1); 
            end
        end
        x(i,1)=x(i,1)/A(i,i);
    end
    x
    %Calculating infinite norm of vector x-xprevious 

    temp=x-xprevious;
    max_temp=temp(1,1);
    for i=2:sizeA
       if abs(temp(i,1))>max_temp
           max_temp=abs(temp(i,1));
       end
    end
end

它实际上适用于 100x100 或更小的矩阵。但是，我的导师希望它适用于 100000x100000 矩阵。起初，甚至很难创建矩阵本身，但我在此处的一些帮助下设法做到了： Matlab Help Center

现在，我以 A 作为参数调用 ex1_3 函数，但它运行起来非常慢。其实它永远不会结束。我怎样才能让它发挥作用？

这是我创建导师想要的特定矩阵的代码： 重要的是它满足以下条件： A(i; i) = 3, A(i - 1; i) = A(i; i + 1) = -1 n=100000

b=ones(100000,1);
b(1,1)=2;
b(100000,1)=2;

i=zeros(299998,1); %Matrix with the lines that we want to put nonzero elements 
j=zeros(299998,1);  %Matrix with the columns that we want to put nonzero elements 
s=zeros(299998,1); %Matrix with the nonzero elements. 
number=1; 
previousNumberJ=0;
numberJ=0;
for k=1:299998 %Our index in i and j matrices
    if mod((k-1),3)==0
        s(k,1)=3;
    else
        s(k,1)=-1;
    end
    if k==1 || k==2
        i(k,1)=1;
        j(k,1)=k;
    elseif k==299997 || k==299998   
        i(k,1)=100000;
        j(k,1)=(k-200000)+2;
    else
        if mod(k,3)==0
            number=number+1;
            numberJ=previousNumberJ+1;
            previousNumberJ=numberJ;
        end
        i(k,1)=number;
        j(k,1)=numberJ;
        numberJ=numberJ+1;
    end
end

A=sparse(i,j,s); %Creating the sparse array

x=ex1_3(A,b);

Answer 1

for 循环在 Matlab 中运行非常缓慢，也许您可​​能想尝试迭代的matrix form：

function x=gseidel(A,b)
    max_temp=10^(-6); %Pass first iteration
    x=b;
    Q=tril(A);
    r=b-A*x;

    for i=1:100
        dx=Q\r; 
        x=x+1*dx; 
        r=b-A*x; 

        % convergence check
       if all(abs(r)<max_temp) && all(abs(dx)<max_temp), return; end
    end

对于您的A 和b，只需 16 步即可收敛。

tril提取A的下三角部分，你也可以在建矩阵的时候得到这个Q。由于Q 已经是三角矩阵，如果不允许使用\ 函数，您可以很容易地solve 等式Q*dx=r。