ヤコビ反復法
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
&ref(ls_jacobi.eq1.gif,nolink,70%);が大きいと，ガウスの消...
以下に示す3元連立方程式を考える．
#ref(ls_jacobi.eq3.gif,nolink,70%)
係数行列の対角成分&ref(ls_jacobi.eq4.gif,nolink,70%);なら...
#ref(ls_jacobi.eq5.gif,nolink,70%)
が導かれる．
&ref(ls_jacobi.eq6.gif,nolink,70%);を初期値として与え，上...
反復回数を&ref(ls_jacobi.eq8.gif,nolink,70%);とし，&ref(l...
n元連立1次方程式の場合は以下となる．
#ref(ls_jacobi.eq10.gif,nolink,70%)
ここで，&ref(ls_jacobi.eq11.gif,nolink,70%);である．
このようにして解を求める方法がヤコビ反復法である．
解が収束したかどうかは&ref(ls_jacobi.eq9.gif,nolink,70%);...
絶対誤差を用いる場合の収束条件は以下となる．
#ref(ls_jacobi.eq13.gif,nolink,70%)
また，相対誤差を用いた場合は，
#ref(ls_jacobi.eq14.gif,nolink,70%)
ここで&ref(ls_jacobi.eq15.gif,nolink,70%);は許容誤差であ...
計算量は反復回数を&ref(ls_jacobi.eq8.gif,nolink,70%);とす...
&ref(ls_jacobi.eq8.gif,nolink,70%);をなるべく小さくするこ...
また，反復法で問題となるのは収束するための条件である．
そもそも収束しないのであれば反復を繰り返しても解が発散す...
収束条件を上記の反復式から導出する．
まず，反復式を以下のような形で表す．
#ref(ls_jacobi.eq17.gif,nolink,70%)
反復回数が十分大きい数&ref(ls_jacobi.eq18.gif,nolink,70%)...
#ref(ls_jacobi.eq19.gif,nolink,70%)
この2つの式から以下が導かれる．
#ref(ls_jacobi.eq20.gif,nolink,70%)
&ref(ls_jacobi.eq21.gif,nolink,70%);は&ref(ls_jacobi.eq22...
&ref(ls_jacobi.eq23.gif,nolink,70%);は&ref(ls_jacobi.eq8....
つまり，反復を繰り返すと1ステップごとに誤差が&ref(ls_jaco...
そのため，解が収束するためには以下の条件を満たす必要があ...
#ref(ls_jacobi.eq25.gif,nolink,70%)
反復式から係数行列の要素を使って表すと，
#ref(ls_jacobi.eq26.gif,nolink,70%)
この式より，係数行列&ref(ls_jacobi.eq27.gif,nolink,70%);...
***ヤコビ反復法の実装 [#bf5910ff]
ヤコビ反復法をC++で実装した例を以下に示す．
#code(C){{
/*!
* ヤコビ反復法(Jacobi iterative method)
* - 解が収束するのは
* ・対角有利(diagonal dominant, 対角要素の絶対値>そ...
* ・係数行列が対称(symmetric)かつ正定(positive defi...
* のどちらかの場合
* @param[inout] A n×nの係数行列とn×1の定数項(b)を併せたn...
* @param[in] n n元連立一次方程式
* @param[inout] max_iter 最大反復数(反復終了後,実際の反...
* @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す)
* @return 1:成功,0:失敗
*/
int JacobiIteration(vector< vector<double> > &A, int n, i...
{
vector<double> x(n, 0.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]); // 絶対誤差の場合
//e += fabs((y[i]-x[i])/y[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;
}
}}
終了行:
&ref(ls_jacobi.eq1.gif,nolink,70%);が大きいと，ガウスの消...
以下に示す3元連立方程式を考える．
#ref(ls_jacobi.eq3.gif,nolink,70%)
係数行列の対角成分&ref(ls_jacobi.eq4.gif,nolink,70%);なら...
#ref(ls_jacobi.eq5.gif,nolink,70%)
が導かれる．
&ref(ls_jacobi.eq6.gif,nolink,70%);を初期値として与え，上...
反復回数を&ref(ls_jacobi.eq8.gif,nolink,70%);とし，&ref(l...
n元連立1次方程式の場合は以下となる．
#ref(ls_jacobi.eq10.gif,nolink,70%)
ここで，&ref(ls_jacobi.eq11.gif,nolink,70%);である．
このようにして解を求める方法がヤコビ反復法である．
解が収束したかどうかは&ref(ls_jacobi.eq9.gif,nolink,70%);...
絶対誤差を用いる場合の収束条件は以下となる．
#ref(ls_jacobi.eq13.gif,nolink,70%)
また，相対誤差を用いた場合は，
#ref(ls_jacobi.eq14.gif,nolink,70%)
ここで&ref(ls_jacobi.eq15.gif,nolink,70%);は許容誤差であ...
計算量は反復回数を&ref(ls_jacobi.eq8.gif,nolink,70%);とす...
&ref(ls_jacobi.eq8.gif,nolink,70%);をなるべく小さくするこ...
また，反復法で問題となるのは収束するための条件である．
そもそも収束しないのであれば反復を繰り返しても解が発散す...
収束条件を上記の反復式から導出する．
まず，反復式を以下のような形で表す．
#ref(ls_jacobi.eq17.gif,nolink,70%)
反復回数が十分大きい数&ref(ls_jacobi.eq18.gif,nolink,70%)...
#ref(ls_jacobi.eq19.gif,nolink,70%)
この2つの式から以下が導かれる．
#ref(ls_jacobi.eq20.gif,nolink,70%)
&ref(ls_jacobi.eq21.gif,nolink,70%);は&ref(ls_jacobi.eq22...
&ref(ls_jacobi.eq23.gif,nolink,70%);は&ref(ls_jacobi.eq8....
つまり，反復を繰り返すと1ステップごとに誤差が&ref(ls_jaco...
そのため，解が収束するためには以下の条件を満たす必要があ...
#ref(ls_jacobi.eq25.gif,nolink,70%)
反復式から係数行列の要素を使って表すと，
#ref(ls_jacobi.eq26.gif,nolink,70%)
この式より，係数行列&ref(ls_jacobi.eq27.gif,nolink,70%);...
***ヤコビ反復法の実装 [#bf5910ff]
ヤコビ反復法をC++で実装した例を以下に示す．
#code(C){{
/*!
* ヤコビ反復法(Jacobi iterative method)
* - 解が収束するのは
* ・対角有利(diagonal dominant, 対角要素の絶対値>そ...
* ・係数行列が対称(symmetric)かつ正定(positive defi...
* のどちらかの場合
* @param[inout] A n×nの係数行列とn×1の定数項(b)を併せたn...
* @param[in] n n元連立一次方程式
* @param[inout] max_iter 最大反復数(反復終了後,実際の反...
* @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す)
* @return 1:成功,0:失敗
*/
int JacobiIteration(vector< vector<double> > &A, int n, i...
{
vector<double> x(n, 0.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]); // 絶対誤差の場合
//e += fabs((y[i]-x[i])/y[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;
}
}}
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
