2012年5月20日日曜日

ガウスの消去法（理論）

ガウスの消去法

ガウスの消去法（ガウスのしょうきょほう）あるいは掃き出し法（はきだしほう）とは、変数の数と方程式の本数が一致した連立一次方程式（線形方程式系）を解く為の解法である。
出典：ガウスの消去法 - Wikipedia(http://ja.wikipedia.org/wiki/ガウスの消去法)

以下，本稿では$n$元線型連立方程式について扱い，$n$本の等式が存在して，解は一意に定まるものとして議論を進める。

１．連立方程式を行列を用いて表す ～拡大係数行列～

線型連立方程式 \[ \left\{\begin{array}{l} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n = b_2 \\ \qquad\vdots \\ a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n = b_n \end{array}\right. \qquad\cdots\cdots\cdots (*) \] を考える。このとき，$n$次正方行列A，$n$次元縦ベクトル$x$，$b$をそれぞれ次のように定義する。 \[ A = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & n_{n2} & \cdots & a_{nn} \end{pmatrix} \quad , \qquad x = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} \quad , \qquad b = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} \] すると，初めの線型連立方程式$(*)$は \[ A x = b \] と表すことができる。
初めの線型連立方程式$(*)$に対して，行列$A$を係数行列といい， \[ (A \mid b) = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & \mid & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & \mid & b_2 \\ \vdots & \vdots & \ddots & \vdots & \mid & \vdots \\ a_{n1} & n_{n2} & \cdots & a_{nn} & \mid & b_n \end{pmatrix}\] によって定義した$(A \mid b)$を拡大係数行列という。
さて，線型連立方程式$(*)$の解が$x_i = c_i ~ (\forall i \in \mathbb{N}, 1 \leq i \leq n)$であるとしよう。
これを言い換えれば，$n$次単位行列$I_n$を用いて， \[ I_n x = c \qquad \left( ただし，c = \begin{pmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{pmatrix} とする。 \right) \] といえる。

２．行列の基本変形と基本行列

ところで，行に関する基本変形を用いて$A \to I_n$と変形することを考えよう。
行に関する基本変形とは以下にあげる3つの変形である。
  • 2つの行を入れ替える
  • ある行を0でない定数倍する
  • ある行に、他のある行の定数倍を加える
$n$次正方行列$P_{ij}$を \[ P_{ij} = \begin{pmatrix} 1 & ~ & ~ & ~ & ~ & ~ & O \\ ~ & \ddots & ~ & ~ & ~ & ~ & ~ \\ ~ & ~ & 0 & ~ & 1 & ~ & ~ \\ ~ & ~ & ~ & \ddots & ~ & ~ & ~ \\ ~ & ~ & 1 & ~ & 0 & ~ & ~ \\ ~ & ~ & ~ & ~ & ~ & \ddots & ~ \\ O & ~ & ~ & ~ & ~ & ~ & 1 \end{pmatrix} \begin{matrix} ~ \\ ~ \\[.5em] \leftarrow i\textrm{行目} \\ ~ \\[.75em] \leftarrow j\textrm{行目} \\ ~ \\ ~ \end{matrix} \] と定義する。これはもちろん$n$次単位行列$I_n$の$i$行目と$j$行目を入れ替えたものに等しい。
ある$n$次正方行列$X$に，左から$P_{ij}$をかけた行列$P_{ij}X$は行列$X$の$i$行目と$j$行目を入れ替えたものに等しくなる。

同様に他の2つの変形についても行列を定義できる。
$n$次正方行列$Q_{i}(c)$を \[ Q_{i}(c) = \begin{pmatrix} 1 & ~ & ~ & ~ & ~ & ~ & ~ & O \\ ~ & \ddots & ~ & ~ & ~ & ~ & ~ \\ ~ & ~ & 1 & ~ & ~ & ~ & ~ \\ ~ & ~ & ~ & c & ~ & ~ & ~ \\ ~ & ~ & ~ & ~ & 1 & ~ & ~ \\ ~ & ~ & ~ & ~ & ~ & \ddots & ~ \\ O & ~ & ~ & ~ & ~ & ~ & 1 \end{pmatrix} \begin{matrix} ~ \\ ~ \\ ~ \\[.5em] \leftarrow i行目 \\ ~ \\ ~ \\ ~ \end{matrix} \] と定義する。これはもちろん$n$単位方行列$I_n$の$i$行$i$列の成分$I(i, i)$を$c$にしたものに等しい。
ある$n$次正方行列$X$に，左から$Q_{i}(c)$をかけた行列$Q_{i}(c)X$は行列$X$の$i$行目を$c$倍したものに等しくなる。

$n$次正方行列$R_{ij}(c)$を \[ R_{ij}(c) = \begin{pmatrix} 1 & ~ & ~ & ~ & ~ & ~ & ~ & O \\ ~ & \ddots & ~ & ~ & ~ & ~ & ~ \\ ~ & ~ & 1 & \cdots & c & ~ & ~ \\ ~ & ~ & ~ & \ddots & \vdots & ~ & ~ \\ ~ & ~ & ~ & ~ & 1 & ~ & ~ \\ ~ & ~ & ~ & ~ & ~ & \ddots & ~ \\ O & ~ & ~ & ~ & ~ & ~ & 1 \end{pmatrix} \begin{matrix} ~ \\ ~ \\[.5em] \leftarrow i\textrm{行目} \\ ~ \\[.75em] \leftarrow j\textrm{行目} \\ ~ \\ ~ \end{matrix} \] と定義する。これはもちろん$n$単位方行列$I_n$の$i$行$j$列の成分$I(i, j)$を$c$にしたものに等しい。
ある$n$次正方行列$X$に，左から$R_{ij}(c)$をかけた行列$R_{ij}(c)X$は行列$X$の$i$行目を$c$倍したものに等しくなる。

以上の用ように$n$次正方行列$X$に左から$P_{ij}$，$Q_{i}(c)$，$R_{ij}(c)$（この3つの行列を総称して基本行列とよぶ。）をかけることで行に関する基本変形を行うことができる。この性質から行に関する基本変形を左基本変形ともいう。
（ちなみに，$n$次正方行列$X$に，右から$P_{ij}$，$Q_{i}(c)$，$R_{ij}(c)$をかけることで列に関する基本変形を行うことができる。この性質から列に関する基本変形を右基本変形ともいう。列に関する基本変形は行に関する基本変形の3つそれぞれについて，行を列と読み替えたものである。）

以上の議論より，当初の目的であった$A \to I_n$の変形は，行に関する基本変形は基本行列を左からかける操作を有限回行うことで達成できる。 （再度述べるが，これはあくまで解が一意に定まるという仮定の下で成り立つ。一般的にはいえない。）

３．ガウスの消去法（掃き出し法）

ガウスの消去法のアルゴリズムは前進消去(forward elimination)と後退代入(backward substitution)と呼ばれる2つのステップよりなる。
以下の疑似コードでは，$a_{i,j}$は$(n ,n+1)$型の拡大係数行列$(A \mid b)$の$i$行$j$列の要素を示す。

３－１．前進消去のアルゴリズム


/* $i$行目の前進消去を行う */
■$Proceduer$ 前進消去($i$: 整数)
｜　　/* 必要に応じてピボット選択 */
｜　　□$If$ $a_{i,i} = 0$ $Then$
｜　　｜　　$j \leftarrow a_{i,p} \neq 0$なる最小のp
｜　　｜　　/* もしそのようなpが存在しなければ，解が一意に定まらないのでエラーとなる */
｜　　｜　　行列$A$の$i$行目と$j$行目を入れ替える
｜　　□$End If$
｜　　$D \leftarrow a_{i,i}$
｜　　/* $i$行目の各要素を$D$で割る */
｜　　□$Loop$ $1 \leq j \leq n$
｜　　｜　　$a_{i,j} \leftarrow \dfrac{a_{i,j}}{D}$
｜　　□$EndLoop$
｜　　/* $i$行目以降の全ての行から$i$行目を引く */
｜　　□$Loop$ $i+1 \leq j \leq n$
｜　　｜　　□$Loop$ $1 \leq k \leq n+1$
｜　　｜　　｜　　$a_{j,k} \leftarrow a_{j,k} - a_{i,k}$
｜　　｜　　□$EndLoop$
｜　　□$EndLoop$
■$EndProceduer$

前進消去では拡大係数行列の下三角の要素を0にし，さらに対角要素を1にすることを行う。
もし，対角要素が行の入れ替えによっても0にしかならないならば，解は一意に定まることができない。


３－２．後退代入のアルゴリズム


/* $i$行目の後退代入を行う */
■$Proceduer$ 後退代入($i$: 整数)
｜　　/* 1ステップ前の解を代入 */
｜　　□$Loop$ $i \geq j \geq 0$
｜　　｜　　$a_{j,n+1} \leftarrow a_{j,n+1} - c_{i+1} \times a_{j,i+1}$
｜　　｜　　$a_{j,i+1} \leftarrow 0$
｜　　□$EndLoop$
｜　　/* $i$行目を用いて1文字消去 */
｜　　□$If$ $a_{i,i} \neq 0$ $Then$
｜　　｜　　$a_{i,n+1} \leftarrow \dfrac{a_{i,n+1}}{a_{i,i}}$
｜　　□$EndIf$
｜　　/* 解を$c_i$とする */
｜　　$c_i \leftarrow a_{i, n+1}$
■$EndProceduer$


３－３．ガウスの消去法全体のアルゴリズム


/* 初期値 */
$a_{i,j} \leftarrow$ 拡大係数行列$(A \mid b)$の$(i, j)$要素
$n \leftarrow $次数
/* 解 */ $c_{i}$ /* 解を入れるベクトル */

/* ガウスの消去法 */
■$Proceduer$ ガウスの消去法
｜　　/* 前進消去 */
｜　　□$Loop$ $1 \leq i \leq n$
｜　　｜　　前進消去$(i)$
｜　　□$EndLoop$
｜　　/* 後退代入 */
｜　　□$Loop$ $n \geq i \geq 1$
｜　　｜　　後退代入$(i)$
｜　　□$EndLoop$
■$EndProceduer$


投稿者 時刻: 0 コメント

2012年5月19日土曜日

ガウスの消去法（実装）

ガウスの消去法の実装を行う。

コードは次のようになる。 
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>

/************************************************************
 * ガウスの消去法による線形方程式の数値計算
 * 
 * (C) 2012- Gofer.
 ************************************************************/

// 方程式の表示
void display_equtation(int N, double *A) {
 int i, j;
 
 for(i=0; i<N; ++i) {
  for(j=0; j<N; ++j) {
   printf("%f x[%d]", A[i*(N+1) + j], j);
   if(j < N-1) printf(" + ");
  }
  printf(" = %f\n", A[i*(N+1) + N]);
 }
 putchar('\n');
 
 return;
}

// n行目とm行目を入れ替える
double* replace_row(int n, int m, int N, double* A) {
 double* line_temp = (double*)malloc(sizeof(double)*(N+1));
 memmove(line_temp,   &A[n*(N+1)], sizeof(double)*(N+1));
 memmove(&A[n*(N+1)], &A[m*(N+1)], sizeof(double)*(N+1));
 memmove(&A[m*(N+1)], line_temp,  sizeof(double)*(N+1));
 free(line_temp);
 return A;
}
 
// i行目の前進消去を行う
double* forward_elimination(int i, int N, double *A) {
 int j, k, J = -1;
 double divide;
 
 // 対角要素が0ならピボット選択を行う
 if(A[i*(N+1) + i] == 0) {
  for(j=i+1; j<N; ++j) {
   if(A[j*(N+1) + i] != 0) {
    J = j;
    break;
   }
  }
  if(J == -1) return NULL;
  replace_row(i, J, N, A);
 }
 
 // 割算する係数を決定
 divide = A[i*(N+1) + i];
 
 // 標準化
 for(j=0; j<(N+1); ++j)
  A[i*(N+1) + j] /= divide;
 
 // i行目以外の行から引く
 for(k=i+1; k<N; ++k) {
  for(j=(N+1)-1; j>=0; --j)
   A[k*(N+1) + j] -= A[k*(N+1) + i] * A[i*(N+1) + j];
 }
 
 return A;
}

// i行目の後退代入
double* backward_substitution(int i, int N, double *A, double *x) {
 // (N-1)-i個の解は決定済み
 int j, k;
 
 // 1ステップ前に求めた解を代入
 k = i+1;
 if(k < N) { 
  for(j=k-1; j>=0; --j) {
   A[j*(N+1) + N] -= x[k] * A[j*(N+1) + k];
   A[j*(N+1) + k] = 0;
  }
 }
 
 // 第i行を用いて1文字消去
 if(A[i*(N+1) + i] != 0) {
  A[i*(N+1) + N] /= A[i*(N+1) + i];
  A[i*(N+1) + i] = 1;
 }
 
 // xベクタに代入
 x[i] = A[i*(N+1) + N];
 
 return A;
}

int main(void) {
 //int N = 2;
 /*
 double A[2][3] = { {1, 2, 5}, 
      {2, 3, 8}};
 */

 //int N = 3;
 /*
 double A[3][4] = { {2, 4, 2, 8}, 
      {4, 10, 3, 17}, 
      {3, 7, 1, 11}};
 */

 //int N = 4;
 /*
 double A[4][5] = { {1, 1, 2, 1, 1}, 
      {1, 1, 3, 2, -2}, 
      {2, -2, 2, -1, 1}, 
      {-1, 1, 0, 1, -1}};
 */
 
 int N = 4;      
 double A[4][5] = { {1, -1, 2, -1, -8}, 
      {2, -2, 3, -3, -20}, 
      {1, 1, 1, 0, -2}, 
      {1, -1, 4, 3, 4}};
 
 int i;
 double* x = (double*)malloc(sizeof(double) * N);
 
 display_equtation(N, (double*)A);
 
 // 前進消去
 for(i=0; i<N; ++i) {
  if(forward_elimination(i, N, (double*)A) == NULL) {
   fprintf(stderr, "Error!\n");
   free(x);
   return -1;
  }
 }
 
 // 後退代入
 for(i=N-1; i>=0; --i)
  backward_substitution(i, N, (double*)A, x);
 
 for(i=0; i<N; ++i)
  printf("x[%d] = %f\n", i, x[i]);
 
 free(x);
 
 return 0;
}
投稿者 時刻: 0 コメント
登録: 投稿 (Atom)