連立一次方程式を解く2

寄稿:東條遼平

LU分解を用いて連立一次方程式を解く

せっかくLU分解を作ったのでそれを用いて連立一次方程式を解いてみます. LU分解を使って連立一次方程式が解けるということはLU分解 のところで書きましたので併せて見てください.例えば係数の行列が3×3であれば,

img174.png

のように下三角行列Lと上三角行列Uに分解でき,

img175.png

という連立一次方程式であれば,

img176.png

を解くことと等価となります.ここでベクトルcを考え,

img177.png

とおくと,

img178.png

となりますが,このベクトルcを求める事を前進代入といいます.cが求まると次に

img177.png

より,ベクトルxを求めることが出来ます.これを後退代入といいます.

ソースコード

行列matrixとベクトルansは引数として受け取った係数行列とyベクトルです.

order = LU(matrix);

for(i=0; i<ans->size; i++){
  if(i != order[i]){
    tempd = ans->v[order[i]];
    ans->v[order[i]] = ans->v[i];
    ans->v[i] = tempd;

    order[order[i]] = order[i];
    order[i] = i;
  }
}

LU分解の計算中に都合がいいように行を入れ換えられるので, どの行とどの行が入れ換えられたかわかるように,例えば0行目と2行目が入れ換えられれば, orderは2,1,0,3,4・・・となっています.そのため,それに合わせるために ベクトルansの方も入れ換えます.

for(i=1; i<height; i++){
  for(j=0; j<i; j++){
    ans->v[i] -= ans->v[j]*matrix->a[i*width + j];
  }
}

それでは前進代入です.下三角行列である行列Lは対角成分が1なので, 他の項を引いてやるだけです.これでベクトルansはベクトルcになりました.

for(i=height-1; i>=0; i--){
  if(!(matrix->a[i*width + i])){
    fprintf(stderr, "Error In FBSubstitution.\n");
    exit(1);
  }

  ans->v[i] /= matrix->a[i*width + i];
  for(j=i-1; j>=0; j--){
    ans->v[j] -= ans->v[i]*matrix->a[j*width + i];
  }
}

ここで後退代入を行っています.後退代入は前進代入と違って対角成分がまちまちなので, 他の項を引いた後に対角成分で割る必要があります.そのためゼロでないかのチェックが必要になります.

Valid XHTML 1.1! home > コンピュータ > プログラミング >
リロード   新規 編集 凍結 差分 添付 複製 改名   トップ 一覧 検索 最終更新 バックアップ   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.5_1 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-jiPowered by PHP 5.3.29HTML convert time to 0.011 sec.