第6回 連立方程式の解法 (1)

6.1 ガウス・ジョルダン(Gauss-Jordan)法

以下のような連立一次方程式(ここでは3元一次を例にとる)
の解を数値演算(解析的ではなく、すなわちコンピュータで)により求める。 ガウス・ジョルダン法は以下の手順により計算を行う。
  1. (1)式を a11 で割り、x1の係数を1にする。→ (1)'式
  2. (2)式から、(1)'式をa21倍したものを引く。→ (2)'式
  3. (3)式から、(1)'式をa31倍したものを引く。→ (3)'式
すると、次のように(2)式、(3)式のx1の係数は0になる。
同様に、
  1. (2)'/a'22 → (2)''
  2. (1)'-(2)''*a'12
  3. (3)'-(2)''*a'32
としてx2の項を消してゆくと、
となる。さらに、
  1. (3)''/a''33 →(3)'''
  2. (1)''-(3)'''*a''13
  3. (2)''-(3)'''*a''23
としてx3の項を消していくと、最終的に
となって解を得ることができる。

6.2 具体例

上の作業を具体的な数値を代入して実行してみる。例えば、次のような3元の連立方程式
の場合、以下のような3行4列の係数行列を作って作業することになる。

6.3 プログラム・コーディング

この操作をプログラムとして実現するには、どうすれば良いだろうか?
また、一般にn元の連立方程式を解く場合には、どのようなプログラムを記述すれば良いだろうか??
プログラミングの際には、以上の操作を 抽象的 あるいは 一般的 に記述する必要があります。

前節の最後の式の左辺が
のように単位行列になっていることに注目しよう。 より一般に、n元の連立方程式を取り扱う場合にも、下図のように単位行列が得られるまで「掃き出し演算」を行えばよいことがわかる。
この作業を一般的に記述すると次のようになる。(一般的な作業手順のことを アルゴリズム と呼ぶ)
  1. ピボットが含まれる行に対し、行の要素(aii, aii+1, ... , ain, bi)をピボット係数(akk)で割る。結果としてピボットは1となる。 なお、ピボットの左側の要素(ak1, ak2, ... , akk-1)は既に0になっているので割らなくて良い。
  2. ピボット行以外の各行に対し、ピボットと同じ列が消去されるように
    (各行の要素)-(ピボット行の要素)×(係数)
    の計算を行う。ここで、係数としてはピボットと同じ列の要素を用いれば、ピボットと同じ列が消去されることになる
  3. 次の行に移る
サンプルプログラム: gauss_jordan.cpp

6.4 例題

  • 上のサンプルプログラムを修正し、以下の4元1次連立方程式の解をガウス・ジョルダン法により求めよ。

    課題 (提出〆切:金曜日22時)

    1. プログラムを改良し、以下の4元1次連立方程式の解をガウス・ジョルダン法により求めよ。
      サンプルプログラムそのままの形では、実行してもプログラムが停止する。 なぜ止まってしまうのかを考え、対策せよ。また、その理由を言葉「ピボット、入れ替え」を用いて説明せよ。
      なおこの説明文は、提出するプログラムファイルの先頭にコメント文として記述すること。 ファイル名は6_1.cpp とすること。
    2. mat06.txt および mat10.txt には、今回の授業で説明した連立方程式を表す係数行列が記述されている。 この txt ファイルをコマンドライン引数から読み込み、連立方程式の解を求めて画面に出力させるプログラムを作成せよ。 ファイル名は 6_2.cpp とすること。
      • ファイルを読み込むまで連立方程式のサイズ(元)がわからない、という条件で実行できるプログラムを作成すること。 すなわち、ファイル読み込みの際に連立方程式のサイズを判定する処理を加える必要がある。
        • 行末やファイル末尾を検出する方法は、検索するといろいろ見つかると思います。
      • 一方で、配列のサイズを動的に調整することはちょっと難しいので(メモリリークしやすい)、 プログラムに含まれる配列のサイズについては、あらかじめ大きめな値で十分な領域を確保しておくこと。 今回は N=20 とする。
        • もちろん、malloc関数やnew関数を使って動的なメモリ管理を行っても構わない。
      • サイズ違いの連立方程式でも正しく動作するかどうかをチェックできるよう、2種類の連立方程式のファイルを提供している。 実行時に読み込むファイル名を変更し、どちらの連立方程式に対しても正しく処理が行われていることを確認すること。
      • 行数と列数の整合性をチェックするエラー処理を含めること。すべての行に対して (列数)=(行数)+1 となっているかがポイントである。

    課題の解答例