ガウスの消去法(Gaussian elimination)と呼ばれる方法では, 前進消去と後退代入という2段階の計算により解を求める. 一般的な(中学校で習う)解き方の一つに加減法というものがある. たとえば,以下のような連立方程式で考える. 加減法では片方の式の両辺に実数をかけて,係数を合わせて,加減算する. ここでは下の式の両辺を2倍して,上の式から引くことにする. あとは,の結果を1番目の式に代入することで,を得る. これを行列で表すと, に対して,2行目にをかけて引くことで, という形にする.そして,一番下の行から未知数を解いて代入していく. 前半の処理では係数行列の対角要素を除いた左下部分がすべて0となる 以下に示すような上三角行列(upper triangle matrix)を求めている. この処理のことを前進消去(forward elimination)と呼ぶ. 前進消去により,未知数はで求められ, このを行目に代入することでが得られ, を行目に代入することでが得られる. これをまで繰り返すことですべての解を得る. この後半部分の処理のことを後退代入(back substitution)と呼ぶ. 前進消去と後退代入により連立1次方程式の解を求める方法をガウスの消去法という. ガウスの消去法の実装†ガウスの消去法を元連立1次方程式に適用する. まず,に定数ベクトルを列の最後に追加したの行列を考える. この行列は拡大行列と呼ばれる. 拡大行列の要素を とする.ここで,C,C++で実装することを考えて要素番号を0始まりにしている. 今後,実装の説明では0で始まるインデックスを用いるので注意. 前進消去は次の式で表される. ここで, 上式によりをで更新することで上三角行列を得る. この処理のC++のコードを以下に示す. // 前進消去(forward elimination) // - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする) for(int k = 0; k < n-1; ++k){ double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが,実際にはj=k+1でよい A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } } 2次元配列A[n][n+1]には係数行列と定数ベクトルが格納されている. 次に後退代入の式を以下に示す. 後退代入のC++のコード例を以下に示す. // 後退代入(back substitution) // - x_nの解はb_n/a_nn,x_nをさらにn-1行の式に代入することでx_(n-1)を求める. // - この作業を1行目まで続けることですべての解を得る. A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; } 計算された結果は別の配列でなくA[i][n]に格納しているので注意. ガウスの消去法の全体のコードは以下. /*! * ガウスの消去法(ピボット交換なし) * @param[inout] A n×nの係数項とn×1の定数項(b)を併せたn×(n+1)の行列.n+1列目に解が入る. * @param[in] n n元連立一次方程式 * @return 1:成功 */ int GaussElimination(vector< vector<double> > &A, int n) { // 前進消去(forward elimination) // - 対角要素をのぞいた左下要素をすべて0にする(上三角行列にする) for(int k = 0; k < n-1; ++k){ double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k; j < n+1; ++j){ // 確認のため左下要素が0になるようにj=kとしたが,実際にはj=k+1でよい A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } } // 後退代入(back substitution) // - x_nの解はb_n/a_nn,x_nをさらにn-1行の式に代入することでx_(n-1)を求める. // - この作業を1行目まで続けることですべての解を得る. A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; } return 1; } 標準の2次元配列の代わりにSTLのvectorを使っている. |