Jacobi法

質問

Jacobi法のアルゴリズム

x[i]_{m+1} = ( b[i] - Σ{j≠i} a[i][j] * x[j]_{m} )
i=1,2,…,n
m=0,1,2,…

を基にプログラム（JAVA）を組んでみました。
（質問に直接関係ないところは省いています）

for( m=0; m<100; m++ ) {

for( i=0; i<n; i++ ) {

for( j=0; j<i; j++) {
if( i != j ) {
tmp = tmp + a[i][j] * x2[j];
}
}
x1[i] = (b[i] - tmp) / a[i][i];
}
}

もう１つ。

for( k=0; k<100; k++ ) {

for( i=0; i<n; i++ ) {
xnew[i] = b[i];
for( j=0; j<n; j++ ) {
if( j != i ) {
xnew[i] = xnew[i] - a[i][j] * xold[j];
}
}
xnew[i] = xnew[i] / a[i][i];
}

2つ目のプログラムの方が正確です。
1つ目のプログラムだと間違いでしょうか？
私には2つ目のプログラムの
xnew[i]=b[i];
の部分がどうしても理解できません。
これを行ってしまうと、アルゴリズムが
x[i]_{m+1} = ( Σ{j≠i}b[i] - a[i][j] * x[j]_{m} )
になってしまわないでしょうか・・・？

どちらでも良いのでしょうか？
納得いきません。

よろしくお願いします。

回答No.1

arrysthmia · ベストアンサー

激しくカテゴリー違いですが…。

(1) xnew[i] = b[i]; について　⇒　それでいいんです。

for( j=0; j<n; j++ ) のループを回る度に、xnew がどのように更新されてゆくか、
追跡してみましょう。

for に入る前　：　xnew[i] = b[i]
j が 0 のとき　：　xnew[i] = b[i] - a[i][0] * xold[0]
j が 1 のとき　：　xnew[i] = b[i] - a[i][0] * xold[0] - a[i][1] * xold[1]
以下同様…

ちゃんと、xnew[i] = b[i] - Σ{j≠i}( a[i][j] * xold[j] )
になっているでしょう？

x[i]_{m+1} = Σ{j≠i}( b[i] - a[i][j] * x[j]_{m} ) だったら、
xnew[i] = b[i]; は、for( j=0; j<n; j++ ) の { } 中に入っているはずです。
for( i が何のループで、for( j が何のループか、確認のこと。

それよりも、a[i][i] で割るところは、xold[i] = xnew[i] / a[i][i]; じゃないとね。
k++ するときに、一回前の xnew を今度は xold と見なすようになるのだから。

(2) １つ目のプログラムについて　⇒　ダメな箇所が三つあります。

・ for( j の範囲が違う。
for( j=0; j<i; j++ ) if( j != i ) では、結局 Σ{j＜i} にしかなりません。
Σ{j≠i}( a[i][j] * x2[j] ) を求めたいなら、２つ目のプログラムのように
for( j=0; j<n; j++ ) if( j != i ) でないと。

・ temp を初期化してない。
for( j=0; j<n; j++ ) の中身を j が 0 で通るとき、未代入の temp を右辺で参照していますね。
これでは、for が終わったとき、temp = (最初に temp に入ってた何か) + Σ{j≠i}( a[i][j] * x2[j] )
になっていまいます。for の前に temp = 0; が必要です。
どうせ初期化するなら、このとき temp = b[i] としてしまえば、２つ目のプログラムの考え方になります。

・ x1 を x2 へコピーしてない。
このままでは、for( m のループを何度回っても、１度だけ回ったのと同じことです。
x1[] と x2[] は、同じ変数でいいんじゃないですかね。

お礼

なるほど、Jacobi法の意味も含めてよく理解できました。
ありがとうございます。