找到雅可比迭代矩阵的谱半径

Question

我正在使用 Matlab 找到 Jacobi 迭代矩阵的光谱半径，其中 A=[4 2 1;1 3 1;1 1 4]。

在 5 次迭代后，我似乎无法输入正确的命令来获取错误的大小。有人可以帮我吗？

以下是我目前在 Matlab 中输入的命令列表：

A=[4 2 1;1 3 1;1 1 4]

A =

 4     2     1
 1     3     1
 1     1     4

 D=diagonal(diagonal(A));L=(A,-1);U=(A,1);

 b=([3 -1 4])

 x0j=zeros([0 0 0]);

 x=D\(-(U+L)*x0j+b);r=b-A*x    %Jacobi iteration.
------------------------------------------------------------------------------
 Error using  * 
 Inputs must be 2-D, o enter code here r at least one input must be scalar.
 To compute element wise TIMES, use TIMES (.*) instead.

Answer 1

矩阵的spectral radius 是其特征值的模的最大值。可以使用max(abs(eig(·))) 简单计算。

但是，正如其他人所注意到的那样，您的整个代码似乎很混乱，实际上并不是有效的 Matlab 代码，所以您的问题并不是真正计算光谱半径，是吗？ algorithm 非常简单易行：

% diagonal part of A and rest
D = diag(diag(A));
R = A - D;

% iteration matrix and offset
T = - inv(D) * R;
C = inv(D) * b;

% spectral radius condition
rho = max(abs(eig(T)));
if rho >= 1
    error('no convergence')
end

% initial guess
x = randn(size(b));

% iteration
while norm(A * x - b) > 1e-15
    x = T * x + C;
end

请注意，我使用inv(D) 直接遵循维基百科中的描述，但可以使用diag(1 ./ diag(D)) 轻松计算对角矩阵的逆矩阵。

我真的不明白为什么需要将R 分成上下两部分。我想这与数值效率有关，但是，Matlab 已经是一种非常有效的矩阵计算高级语言。所以实际上，当你可以简单地写 A \ b 时，实际上没有必要在其中显式地实现 Jacobi 算法——我猜是为了教育目的。