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을 사용할 때 언어의 기능을 사용해보십시오. 귀하가 여전히 낮은 수준의 언어를 사용하고있는 것처럼 가장하지 마십시오. 따라서, 우리는

D가의 대각선을 포함하는 대각 행렬이다

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

로서 코비 반복을 쓸 수 있으며, R은 오프 - 대각선 엘리먼트들의 행렬이므로, 대각선에 제로가있다. 우리는 어떻게 MATLAB에서이 작업을 수행 할 수 있습니다?

먼저 D와 R을 간단한 방법으로 작성하십시오.

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

이제 대각선 행렬의 역행렬 계산은 어리 석다는 것을 알아야합니다. 더 대각선의 각 엘리먼트의 역을 계산하는 것이다.

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

그래서, 지금 우리가

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

같은 단일 코비가 반복 쓸 수는 NO 중첩 루프가 필요되지 않았 음을 참조하십시오. 우리는 이들 행렬에 색인을 귀찮게하지 않습니다. 이제 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에 대한 더 나은 시작 값을 선택했을 것입니다. 그리고 마지막으로,보다 완벽한 도움말과 참조를 제공하고자합니다.

보고 싶은 것은 좋은 코드의 일반적인 특성입니다.