1

疑似コードを MatLab アルゴリズムに変換する際にいくつかの問題があります。具体的には、やり方がよくわからない部分が 1 つあります。不確かなところまで疑似コードを書きます。

 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)含まれないようにするにはどうすればよいですか?

どんな助けでも大歓迎です！

4

2 に答える 2

4

あなたはコードを提供したので、こうして努力しました...私はより良いアプローチを指摘します。MATLABを使用するときは、その言語の機能を使用してみてください。まだ低水準言語を使用しているふりをしないでください。したがって、ヤコビ反復を次のように書くことができます。

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));

これで、次のように1つのJacobi反復を記述できます。

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のより適切な開始値を選択した可能性があります。最後に、参照を使用して、より完全なヘルプを提供したいと思います。

あなたが見たいのは、良いコードの一般的な特徴です。

于 2012-09-11T21:00:43.810 に答える
2

明示的なインデックス配列を作成し、不要なインデックスを削除できます

% inside of the loop
idx = 1:n;
idx(i) = [];
u(i) = (b(i)-sum(A(i,idx).*x(idx)))/A(i,i);
于 2012-09-11T20:10:28.397 に答える