3 連立一次方程式

3.1 ガウス・ジョルダン法の基本的な考え方

ガウス・ジョルダン(Gauss-Jordan)法というのは、連立方程式 (6)を次にように変形させて、解く方法である。

$$\begin{aligned}\begin{pmatrix}1 & 0 & 0 & \cdots & 0 \\ 0 & 1 &... ... \\ b_3^\prime \\ \vdots\\ b_N^\prime \end{pmatrix} \end{aligned}$$

この式から明らかに、求める解 $x_{i}=b_{i}^\prime$となる。これをどう やって求めるか?。コンピューターで実際に計算する前に、人力でガウス・ジョ ルダン法で計算してみる。例として、

$$\left\{ \begin{aligned}3x+2y+z&=10\\ x+y+z&=6\\ x+2y+2z&=11 \end{aligned} \right.$$

をガウス・ジョルダン法で解を求める。

まずは、1行目の$x$の係数を1に、2と3行目のそれは0にします。そのために、 1行目は$x$の係数の値で割る。2行目と3行目は、1行目に適当な係数を掛けて引く。次の ようにする。

$$\left\{ \begin{aligned}x+\frac{2}{3}y+\frac{1}{3}z&=\frac{10}{3... ...\ 0+\frac{4}{3}y+\frac{5}{3}z&=\frac{23}{3} \end{aligned} \right.$$

つぎに、2行目の$y$の係数を1にして、1と3行目のそれを0にする。

$$\left\{ \begin{aligned}x+\frac{2}{3}y+\frac{1}{3}z&=\frac{10}{3... ...n{aligned}x+0-z&=-2\\ 0+y+2z&=8\\ 0+0-z&=-3 \end{aligned} \right.$$

同じことを$z$についても繰り返えす。すると、

$$\left\{ \begin{aligned}x+0-z&=-2\\ 0+y+2z&=8\\ 0+0+z&=3 \end{al... ...egin{aligned}x+0+0&=1\\ 0+y+0&=2\\ 0+0+z&=3 \end{aligned} \right.$$

となる。従って、連立方程式(24)の解は、

$$\left\{ \begin{aligned}x=1 \\ y=2 \\ z=3 \end{aligned} \right.$$

となる。これがガウス・ジョルダン法である。もっともらしい名前が付けられているが、 大したことなはい。

これを係数行列とベクトルで書くと、次のようになる。

3.2 ガウス・ジョルダン法のC言語の関数

ピボット選択は行わないで、逆行列も求めないのガウス・ジョルダン法で連立方程式を計 算するプログラムを示す。このプログラムの動作は、次の通りである。

• 仮引数「n」は、解くべき連立方程式の未知数の数である。
• 仮引数の配列「a」と「b」は、係数行列 $\boldsymbol{A}$と非同次項 $\boldsymbol{b}$である。 値は、呼び出し元からにアドレス渡しで送られる。
  • 係数行列は、配列「a[1][1]」〜「a[n][n]」に格納されている。
  • 非同次項は、配列「b[1]」〜「b[n]」に格納されている。
• 連立方程式の解 $\boldsymbol{x}$は、配列「b[1]」〜「b[n]」に格納される。
• このプログラムでの処理が終了すると、配列「a[1][1]」〜「a[n][n]」は単位 行列になる。。

	/* ========== ガウスジョルダン法の関数 =================*/
	void gauss_jordan(int n, double a[][100], double b[]){

	   int ipv, i, j;
	   double inv_pivot, temp;

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 対角成分=1(ピボット行の処理) ---- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }
	}

ホームページ: Yamamoto's laboratory
著者: 山本昌志