matlab - 根据输入在matlab中实现的异常行为算法

Question

我正在做科学计算的家庭作业，特别是 matlab 中的迭代方法 Gauss-Seidel 和 SOR，问题是对于一个矩阵给我意想不到的结果(解决方案不收敛)而另一个矩阵收敛。

这里是sor的代码，其中:

• A:系统的矩阵A * x = b
• Xini:初始迭代数组
• b: 与系统无关的数组 A * x = b
• maxiter:最大迭代次数
• tol:公差；
• 特别是，SOR 方法将接收称为 w 的第六个参数，它对应于松弛参数。

这是 sor 方法的代码:

 function [x2,iter] = sor(A,xIni, b, maxIter, tol,w)            
     x1 = xIni;
     x2 = x1;
     iter = 0;
     i = 0;
     j = 0;
     n = size(A, 1);

     for iter = 1:maxIter,
         for i = 1:n
             a = w / A(i,i);
             x = 0;
             for j = 1:i-1
                 x = x + (A(i,j) * x2(j));
             end
             for j = i+1:n
                 x = x + (A(i,j) * x1(j));
             end
             x2(i) = (a * (b(i) - x)) + ((1 - w) * x1(i)); 
         end
         x1 = x2;
         if (norm(b - A * x2) < tol);
             break;
         end
     end

这是高斯赛德尔法的代码:

function [x, iter] = Gauss(A, xIni, b, maxIter, tol)

x = xIni;
xnew = x;
iter = 0;
i = 0;
j = 0;
n = size(A,1);

for iter = 1:maxIter,
    for i = 1:n
        a = 1 / A(i,i);
        x1 = 0;
        x2 = 0;
        for j = 1:i-1
            x1 = x1 + (A(i,j) * xnew(j));
        end
        for j = i+1:n
            x2 = x2 + (A(i,j) * x(j));
        end
        xnew(i) = a * (b(i) - x1 - x2);
    end
    x= xnew;
    if ((norm(A*xnew-b)) <= tol);
        break;
    end
end

对于这个输入:

A = [1 2 -2; 1 1 1; 2 2 1];
b = [1; 2; 5];

当调用函数 Gauss-Seidel 或 sor 时:

[x, iter] = gauss(A, [0; 0; 0], b, 1000, eps)
[x, iter] = sor(A, [0; 0; 0], b, 1000, eps, 1.5)

高斯的输出是:

x =

  1.0e+304 *

    1.6024
   -1.6030
    0.0011


 iter =

            1000

对于 sor 来说是:

x =

   NaN
   NaN
   NaN


iter =

        1000

但是对于以下系统能够找到解决方案:

A = [    4 -1  0 -1  0  0;
        -1  4 -1  0 -1  0;
         0 -1  4  0  0 -1;
        -1  0  0  4 -1  0;
         0 -1  0 -1  4 -1;
         0  0 -1  0 -1  4   ]
b = [1 0 0 0 0 0]'

解决方案:

[x, iter] = sor(A, [0; 0; 0], b, 1000, eps, 1.5)
x =

    0.2948
    0.0932
    0.0282
    0.0861
    0.0497
    0.0195


iter =

    52

方法的行为取决于两个矩阵的条件？因为我注意到第二个矩阵比第一个矩阵条件更好。有什么建议吗？

Answer 1

来自wiki article on Gauss-Seidel :

convergence is only guaranteed if the matrix is either diagonally dominant, or symmetric and positive definite

由于 SOR 与 Gauss-Seidel 相似，我希望 SOR 也适用相同的条件，但您可能需要查一下。

您的第一个矩阵绝对不是对角占优或对称的。但是，您的第二个矩阵是对称且正定的(因为 all(A==A.') 和 all(eig(A)>0))。

如果您使用 Matlab 的默认方法 (A\b) 作为“真实”解，并绘制每次迭代与“真实”解之间差异的范数，那么您将得到下面两张图。很明显，第一个矩阵永远不会收敛，而第二个矩阵在几次迭代后已经产生了可接受的结果。

在野外应用它们之前始终了解您的算法的局限性:)