3

我正在尝试在 MATLAB 中实现 Jacobi 迭代，但无法使其收敛。我已经在网上和其他地方查看了用于比较的工作代码，但我找不到任何与我的代码相似并且仍然有效的代码。这是我所拥有的：

function x = Jacobi(A,b,tol,maxiter)

n = size(A,1);
xp = zeros(n,1); 
x = zeros(n,1); 
k=0; % number of steps

while(k<=maxiter)
    k=k+1;

    for i=1:n
       xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:n)*x(i+1:n));
    end

    err = norm(A*xp-b);

    if(err<tol)
        x=xp;
        break;
    end

    x=xp;

end

无论我使用什么 A 和 b，这都会爆炸。这可能是我忽略的一个小错误，但如果有人能解释什么是错误的，我将非常感激，因为这应该是正确的，但在实践中并非如此。

4

3 回答 3

7

你的代码是正确的。它似乎不起作用的原因是因为您指定的系统在使用 Jacobi 迭代时 可能不会收敛。

具体来说（感谢@Saraubh），如果您的矩阵严格对角占优，则此A方法将 收敛。换句话说，对于矩阵中的每一行，没有对角线系数的行i中所有列的绝对总和必须小于对角线本身。换句话说：jii

废话

但是，即使不满足此条件，也有一些系统会与 Jacobi 收敛，但在尝试将 Jacobi 用于您的系统之前，您应该将此作为一般规则。如果你使用 Gauss-Seidel，它实际上更稳定。唯一的区别是，当您在行中前进时，您正在重新使用 的解决方案x并将其输入到其他变量中。要制作这个 Gauss-Seidel，您所要做的就是在for循环中更改一个字符。从这里改变它：

xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:n)*x(i+1:n));

对此：

xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*xp(1:i-1) - A(i,i+1:n)*x(i+1:n));
                                    **HERE**

以下是我将向您展示的两个示例：

  1. 我们指定了一个不通过 Jacobi 收敛的系统，但有一个解决方案。这个系统不是对角占优的。
  2. 我们指定一个确实由 Jacobi 收敛的系统。具体来说，这个系统是对角线占优的。

示例 #1

A = [1 2 2 3; -1 4 2 7; 3 1 6 0; 1 0 3 4];
b = [0;1;-1;2];
x = Jacobi(A, b, 0.001, 40)
xtrue = A \ b

x =

1.0e+09 *

4.1567
0.8382
1.2380
1.0983

xtrue =

-0.1979
-0.7187
 0.0521
 0.5104

现在，如果我使用 Gauss-Seidel 解决方案，这就是我得到的：

x =

-0.1988
-0.7190
0.0526
0.5103

哇！它收敛于 Gauss-Seidel 而不是 Jacobi，即使该系统不是对角占优的，我可能对此有一个解释，稍后我会提供。

示例 #2

A = [10 -1 2 0; -1 -11 -1 3; 2 -1 10 -1; 0 3 -1 8];
b = [6;25;-11;15];
x = Jacobi(A, b, 0.001, 40);
xtrue = A \ b

x =

0.6729
-1.5936
-1.1612
2.3275

xtrue =

0.6729
-1.5936
-1.1612
2.3274

这就是我从 Gauss-Seidel 得到的：

x =

0.6729
-1.5936
-1.1612
2.3274

这对于两者来说肯定是收敛的，并且系统是对角线占主导地位的。

因此，您的代码没有任何问题。您只是指定了一个无法使用 Jacobi 解决的系统。对于围绕这种求解的迭代方法，最好使用 Gauss-Seidel。原因是因为您立即使用来自当前迭代的信息并将其传播到其余变量。Jacobi 没有这样做，这就是它更快发散的原因。对于 Jacobi，您可以看到示例 #1 未能收敛，而示例 #2 收敛。Gauss-Seidel 对两者都收敛了。事实上，当它们都收敛时，它们非常接近真正的解决方案。

同样，您需要确保您的系统在对角线上占主导地位，以保证您具有收敛性。不执行此规则...嗯...您将承担风险，因为它可能会或可能不会收敛。

祝你好运！

于 2014-07-14T09:07:31.310 回答
1

虽然这并没有指出您的代码中的问题，但我相信您正在寻找Numerical Methods: Jacobi File Exchange Submission。

%JACOBI   Jacobi iteration for solving a linear system.
% Sample call
%   [X,dX] = jacobi(A,B,P,delta,max1)
%   [X,dX,Z] = jacobi(A,B,P,delta,max1)

它似乎完全按照您的描述进行。

于 2014-07-14T08:30:35.653 回答
0

正如其他人指出的那样，并非所有系统都使用 Jacobi 方法收敛，但他们没有指出为什么？实际上只有一小部分系统与 Jacobi 方法收敛。

收敛标准是“一行中所有系数（非对角线）的总和”必须小于“该行对角线位置的系数”。所有行都必须满足此条件。您可以在以下位置阅读更多信息：Jacobi Method Convergence

在决定使用 Jacobi 方法之前，您必须查看数值方法是否满足此标准。Gauss-Seidel 方法具有稍微宽松的收敛标准，允许您将其用于大多数有限差分类型的数值方法。

于 2014-07-15T06:51:33.140 回答