Jacobi方法的MatLab算法

Question

我在將僞代碼轉換爲MatLab算法時遇到了一些困難。具體來說，我不知道該怎麼做。我會寫的僞代碼最多的地方我不清楚這一點：

input n, (a_{ij}), (b_i), (x_i), M 
    for k = 1 to M do 
     for i = 1 to n do 

      u_i = (b_i - sum[(from j = 1, (j ≠ i), to n) a_{ij} * x_j])/a_{ii} 
     end do 

我有難度的是，當我不得不寫的總和部分。我想不出的是如何編寫算法，使得j≠i的術語不包含在內。到目前爲止，我已經寫道：

function [k,x] = jacobimethod(A,b,M) 
n = length(A); 
u = zeros(1,n); 
x = zeros(1,n); 
for k = 1:M 
    for i = 1:n 
     u(i) = (b(i) - (A(i 

而這就是我卡住的地方。到目前爲止，我所有的算法都涉及應該包含所有項的總和，甚至包括j = i的項。在這種情況下，加和項將只是：

A(i,1:n)*x(1:n)' 

但我怎麼能修改此使A(i,i)不包括？

任何幫助將不勝感激！

Answer 1

你可以做一個明確的指標陣列和刪除索引你不需要

% inside of the loop 
idx = 1:n; 
idx(i) = []; 
u(i) = (b(i)-sum(A(i,idx).*x(idx)))/A(i,i); 

Answer 2

既然你提供的代碼，從而作出了努力。我會指出一個更好的方法。當您使用MATLAB時，請嘗試使用該語言的功能。不要假裝你還在使用低級語言。因此，我們可以寫一個Jacobi迭代作爲

X_(n+1) = inv(D)*(b - R*X_n) 

其中d是包含對角線A的對角線矩陣，R是A的非對角元素的矩陣，所以有在對角線上的零。我們如何在MATLAB中做到這一點？

首先，以簡單的方式構建D和R.

D = diag(diag(A)); 
R = A - D; 

現在，我們應該認識到，計算對角矩陣的逆是愚蠢的。更好的辦法是計算對角線上每個元素的倒數。

Dinv = diag(1./diag(A)); 

所以，現在我們可以寫一個單一的雅可比迭代，因爲

X = Dinv*(b - R*X); 

見人說人不需要嵌套循環。我們根本不打算索引這些矩陣。現在將它們全部包裝在一個MATLAB函數中。保持友善，檢查問題，並大膽使用評論。

============================================== ====

function [X,residuals,iter] = JacobiMethod(A,b) 
% solves a linear system using Jacobi iterations 
% 
% The presumption is that A is nxn square and has no zeros on the diagonal 
% also, A must be diagonally dominant for convergence. 
% b must be an nx1 vector. 

n = size(A,1); 
if n ~= size(A,2) 
    error('JACOBIMETHOD:Anotsquare','A must be n by n square matrix') 
end 
if ~isequal(size(b),[n 1]) 
    error('JACOBIMETHOD:incompatibleshapes','b must be an n by 1 vector') 
end 

% get the diagonal elements 
D = diag(A); 
% test that none are zero 
if any(D) == 0 
    error('JACOBIMETHOD:zerodiagonal', ... 
    'The sky is falling! There are zeros on the diagonal') 
end 

% since none are zero, we can compute the inverse of D. 
% even better is to make Dinv a sparse diagonal matrix, 
% for better speed in the multiplies. 
Dinv = sparse(diag(1./D)); 
R = A - diag(D); 

% starting values. I'm not being very creative here, but 
% using the vector b as a starting value seems reasonable. 
X = b; 
err = inf; 
tol = 100*eps(norm(b)); 
iter = 0; % count iterations 
while err > tol 
    iter = iter + 1; 
    Xold = X; 

    % the guts of the iteration 
    X = Dinv*(b - R*X); 

    % this is not really an error, but a change per iteration. 
    % when that is stable, we choose to terminate. 
    err = norm(X - Xold); 
end 

% compute residuals 
residuals = b - A*X; 

======================================= ===========

讓我們看看它是如何工作的。

A = rand(5) + 4*eye(5); 
b = rand(5,1); 
[X,res,iter] = JacobiMethod(A,b) 

X = 
     0.12869 
    -0.0021942 
     0.10779 
     0.11791 
     0.11785 

res = 
    5.7732e-15 
    1.6653e-14 
    1.5654e-14 
    1.6542e-14 
    1.843e-14 

iter = 
    39 

它是否收斂到我們從反斜槓獲得的解決方案？

A\b 
ans = 
     0.12869 
    -0.0021942 
     0.10779 
     0.11791 
     0.11785 

這對我來說很好。更好的代碼可能會檢查對角線優勢以試圖預測代碼何時失敗。我可能會選擇對解決方案更加智能的容忍度，或者對X更好的起始值。最後，我想提供更完整的幫助和參考。

你想看到的是良好代碼的一般特性。