Jacobi求解器进入无限循环

Question

我似乎找不到解决无限循环的办法。我编写了一个Jacobi求解器来求解线性方程组。

这是我的代码：

function [x, i] = Jacobi(A, b, x0, TOL)
[m n] = size(A);
i = 0;
x = [0;0;0];
while (true)
    i =1;

     for r=1:m
        sum = 0;
        for c=1:n
            if r~=c
                sum = sum + A(r,c)*x(c);
            else
                x(r) = (-sum + b(r))/A(r,c);
            end    
            x(r) = (-sum + b(r))/A(r,c);
 xxx   end                                       xxx
    end
    if abs(norm(x) - norm(x0)) < TOL;
        break
    end
    x0 = x;
    i = i + 1;
end

当我终止代码时，它以xxx行结尾。

Answer 1

代码不能工作的原因是if语句在for循环中的逻辑。具体来说，您需要先积累不属于该行对角线的特定行的所有值。一旦完成，您就可以然后执行除法了。您还需要确保将您正在关注的行的对角线系数除以A的对角线系数，该系数对应于您正试图求解的x的组件。您还需要在循环开始时删除i=1语句。你每次都要重置i。

换言之：

function [x, i] = Jacobi(A, b, x0, TOL)
[m n] = size(A);
i = 0;
x = [0;0;0];
while (true)
     for r=1:m
        sum = 0;
        for c=1:n
            if r==c %// NEW
                continue;
             end
            sum = sum + A(r,c)*x(c); %// NEW
        end                                      
       x(r) = (-sum + b(r))/A(r,r); %// CHANGE
    end
    if abs(norm(x) - norm(x0)) < TOL;
        break
    end
    x0 = x;
    i = i + 1;
end

示例使用：

A = [6 1 1; 1 5 3; 0 2 4]
b = [1 2 3].';
[x,i] = Jacobi(A, b, [0;0;0], 1e-10)

x =

   0.048780487792648
  -0.085365853612062
   0.792682926806031


i =

    20

这意味着需要20次迭代才能实现具有容差1e-10的解决方案。将此与MATLAB的内置逆进行比较：

x2 = A \ b

x2 =

   0.048780487804878
  -0.085365853658537
   0.792682926829268

正如您所看到的，我指定了对1e-10的容忍度，这意味着我们保证有10位小数点的精度。我们当然可以看到10位小数点之间的精度雅各比给我们与MATLAB给我们的内置。