[線代] Jacobi iteration

[simon860730 (近物GG)]

這是一題我們數據分析的作業
要用Jacobi iteration method解方程式(?)

要解的東西如下

x1 +  x2 +   x3 +   x4 +   x5 =  5
x1 + 2x2 +  3x3 +  4x4 +  5x5 =  3
x1 + 3x2 +  6x3 + 10x4 + 15x5 =  1
x1 + 4x2 + 10x3 + 20x4 + 35x5 = -3
x1 + 5x2 + 15x3 + 35x4 + 70x5 = -5

但是我怎麼做怎麼發散...
還請大家幫我看一下
是這個用Jacobi下去做原本就會發散
(有辦法轉換矩陣的型態使其不發散嗎?)
還是我的code打得有問題

我是寫MATLAB：

A = [1 1  1  1  1 ;
     1 2  3  4  5 ;
     1 3  6 10 15 ;
     1 4 10 20 35 ;
     1 5 15 35 70 ];

b = [5;3;1;-3;-5];

n = sqrt(numel(A));
D = zeros(n);
R = zeros(n);
x = zeros(n,1);

% 把A拆解成D跟R

for i = 1:n
    for j = 1:n
        if j == i
            D(i,j) = A(i,j);
        else
            R(i,j) = A(i,j);
        end
    end
end

itr = 0;

while itr < 100
    xnext = inv(D)*(b-R*x);
    x = xnext
    itr =itr + 1
end

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.115.204.20
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1512889886.A.ED9.html
推 LiamIssac   : maybe you should try stopping the algorithm as x
→ LiamIssac   :  and xnext have very small deviation (not using
→ LiamIssac   :  100 iterations)
when it iterate 100 times , I get my answer
x =

   1.0e+38 *

   -4.7250
   -5.4609
   -3.5379
   -1.8584
   -0.8646

when I try to calculate the answer with Gauss-Seidel iteration

x =

   14.9999
  -31.9998
   41.9997
  -25.9998
    6.0000

itr = 4295
※ 編輯: simon860730 (140.115.204.20), 12/10/2017 15:26:01
推 LPH66       : 我看維基百科是說 Jacobi method 要求 A 是對角優勢
→ LPH66       : 也就是主對角線上的絕對值要在橫列中最大
→ LPH66       : 但這個方程組看起來似乎很難轉成這種矩陣...
推 sin55688    : 先檢查跌代矩陣的eig是否滿足收斂條件
推 sin55688    : 就是 inv(D)R 譜半徑要小於1 才會收斂
inv(D)*R =

         0   -1.0000   -1.0000   -1.0000   -1.0000
   -0.5000         0   -1.5000   -2.0000   -2.5000
   -0.1667   -0.5000         0   -1.6667   -2.5000
   -0.0500   -0.2000   -0.5000         0   -1.7500
   -0.0143   -0.0714   -0.2143   -0.5000         0

eig(inv(D)*R) =

   -2.4808
   -0.1831
    0.6999
    0.9654
    0.9986

爆炸了  不會收斂呢
※ 編輯: simon860730 (140.115.204.20), 12/10/2017 16:27:02
※ 編輯: simon860730 (140.115.204.20), 12/10/2017 16:48:44
※ 編輯: simon860730 (140.115.204.20), 12/10/2017 16:50:15