ヤコビ反復法†
が大きいと，ガウスの消去法では時間がかかりすぎる(
)．これに対して，ガウスの消去法のように直接解くのではなく，適当な初期値から出発して計算を反復することで解に近づけていく方法を反復法と呼ぶ．ここでは最も基本的な反復法であるヤコビ(Jacobi)反復法について述べる．
以下に示す3元連立方程式を考える．
係数行列の対角成分
ならば，
が導かれる．
を初期値として与え，上式に基づき
と計算していく．
反復回数を
とし，を十分大きくすると
は解に収束する．
n元連立1次方程式の場合は以下となる．
ここで，
である．
このようにして解を求める方法がヤコビ反復法である．
解が収束したかどうかはと
の値を比べて，その差が許容誤差内に収まったかどうかで判断する．
絶対誤差を用いる場合の収束条件は以下となる．
また，相対誤差を用いた場合は，
ここで
は許容誤差であり，ユーザが任意の値を設定する．
計算量は反復回数をとすると
となり，
をなるべく小さくすることが計算量を減らすために重要となる．
また，反復法で問題となるのは収束するための条件である．
そもそも収束しないのであれば反復を繰り返しても解が発散するだけである．
収束条件を上記の反復式から導出する．
まず，反復式を以下のような形で表す．
反復回数が十分大きい数
であるとき，以下の式が成り立つ．
この2つの式から以下が導かれる．
は
ステップでの誤差，
はステップでの誤差である．
つまり，反復を繰り返すと1ステップごとに誤差が
倍になるということである．
そのため，解が収束するためには以下の条件を満たす必要がある．
反復式から係数行列の要素を使って表すと，
この式より，係数行列
の各行において対角要素の絶対値が非対角要素の絶対値の和より大きい場合に収束する．このような条件を対角優位(diagonal dominant)と呼ぶ．
ヤコビ反復法の実装†
ヤコビ反復法をC++で実装した例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
| |
int JacobiIteration(vector< vector<double> > &A, int n, int &max_iter, double &eps)
{
vector<double> x(n, 0.0); vector<double> y(n, 0.0);
double e = 0.0;
int k;
for(k = 0; k < max_iter; ++k){
for(int i = 0; i < n; ++i){
y[i] = A[i][n];
for(int j = 0; j < n; ++j){
y[i] -= (j != i ? A[i][j]*x[j] : 0.0);
}
y[i] /= A[i][i];
}
e = 0.0;
for(int i = 0; i < n; ++i){
e += fabs(y[i]-x[i]); }
if(e <= eps){
break;
}
swap(x, y);
}
max_iter = k;
eps = e;
for(int i = 0; i < n; ++i){
A[i][n] = y[i];
}
return 1;
}
|
ガウス・ザイデル反復法†
ヤコビ反復法では
ステップでの
だけを使って
での値を求めていた．
しかし，
と順番に計算した場合，
を求めるときには
が，
を求めるときには
がすでに計算されている．
これら最新の値を使うのがガウス・ザイデル反復法(Gauss-Seidel iteration method)である．
たとえば，3元連立1次方程式の場合，反復式は以下となる．
収束判定条件はヤコビ法と同じであり，
収束するための係数行列の条件も同じである．
ただし，収束までの反復回数はヤコビ法よりも少なくすむ．
ガウス・ザイデル反復法の実装†
ヤコビ反復法では
と
の2つを格納するために2つの配列を用意したが，
ガウス・ザイデル法では常に最新の値を用いるため配列は1つでよい．
ガウス・ザイデル反復法をC++で実装した例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
| |
int GaussSeidel(vector< vector<double> > &A, int n, int &max_iter, double &eps)
{
vector<double> x(n, 0.0); double tmp;
double e = 0.0;
int k;
for(k = 0; k < max_iter; ++k){
e = 0.0;
for(int i = 0; i < n; ++i){
tmp = x[i];
x[i] = A[i][n];
for(int j = 0; j < n; ++j){
x[i] -= (j != i ? A[i][j]*x[j] : 0.0);
}
x[i] /= A[i][i];
e += fabs(tmp-x[i]); }
if(e <= eps){
break;
}
}
max_iter = k;
eps = e;
for(int i = 0; i < n; ++i){
A[i][n] = x[i];
}
return 1;
}
|
SOR法†
ヤコビ法やガウス・ザイデル法において各ステップで計算された
と
を用いて，さらに収束を加速させる方法が逐次加速緩和法(SOR法:Successive Over-Relaxation method)である．
各ステップにおいて補正量
を用い，次の式で次ステップの値
を計算する．
ここで，
は加速係数であり，通常，
である．
SOR法を使うとの値によっては解の収束を速めることができる．
ただし，問題にも依存するため最適な加速係数の算出が難しい．
SOR法の実装†
SOR法をC++で実装した例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
| |
int SOR(vector< vector<double> > &A, int n, double w, int &max_iter, double &eps)
{
vector<double> x(n, 0.0); double tmp;
double e = 0.0;
int k;
for(k = 0; k < max_iter; ++k){
e = 0.0;
for(int i = 0; i < n; ++i){
tmp = x[i];
x[i] = A[i][n];
for(int j = 0; j < n; ++j){
x[i] -= (j != i ? A[i][j]*x[j] : 0.0);
}
x[i] /= A[i][i];
x[i] = tmp+w*(x[i]-tmp);
e += fabs(tmp-x[i]); }
if(e <= eps){
break;
}
}
max_iter = k;
eps = e;
for(int i = 0; i < n; ++i){
A[i][n] = x[i];
}
return 1;
}
|
参考文献†
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
