作業中のメモ

よく「計算機」を使って作業をする．知らなかったことを中心にまとめるつもり．

トップ > 数学 > ガウスの消去法のプログラム

ガウスの消去法のプログラム

数学 プログラミング

前回は，関数近似を行うことを考えた．このとき，係数を求めるために，連立方程式を解く必要が生じた． 今回は，とりあえず，ガウスの消去法を用いて，連立方程式を解くことを考える．

ガウスの消去法

ガウスの消去法は，中学，高校のときに習った連立方程式の 1 つである消去法を用いる． ガウスの消去法では，行列表示された連立方程式の係数行列を上三角行列に変形（前進消去）し，今度は逆方向に解を求めていく（後退代入）という手法をとる．

例題

例えば，以下のような 3 × 3 の係数行列と右辺ベクトルから解を求めることを考える．

 { \displaystyle
  \begin{equation}
    \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        1 & 1 & -1\\
        1 & -2 & 3\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        4\\
        -1\\
      \end{array}
    \right)
  \end{equation}
}

この行列を上三角行列に変形したいので，まずは，1 列目に関して， {\dfrac{(2, 1)\mbox{成分}}{(1, 1)\mbox{成分}}} の値を計算し， 2 行目から {\dfrac{(2, 1)\mbox{成分}}{(1, 1)\mbox{成分}} \times \mbox{1 行目}} を引けばよい．そして，右辺ベクトルにも同様のことを行う．この場合， {\dfrac{(2, 1)\mbox{成分}}{(1, 1)\mbox{成分}} = \dfrac{1}{2}} であるので，

 { \displaystyle
  \begin{eqnarray}
    & & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        1 & 1 & -1\\
        1 & -2 & 3\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        4\\
        -1\\
      \end{array}
    \right)\\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        1 - \dfrac{1}{2} \cdot 2 & 1 - \dfrac{1}{2} \cdot 1 & -1 - \dfrac{1}{2} \cdot (-2)\\
        1 & -2 & 3\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        4 - \dfrac{1}{2} \cdot 1\\
        -1\\
      \end{array}
    \right) \\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        1 & -2 & 3\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        -1\\
      \end{array}
    \right) 
  \end{eqnarray}
}

となる．同様に， {(3, 1)} 成分に対しても計算を行う．この場合， {\dfrac{(3, 1)\mbox{成分}}{(1, 1)\mbox{成分}} = \dfrac{1}{2}} であるので，

 { \displaystyle
  \begin{eqnarray}
    & & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        1 & -2 & 3\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        -1\\
      \end{array}
    \right) \\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        1 - \dfrac{1}{2} \cdot 2 & -2 - \dfrac{1}{2} \cdot 1 & 3 - \dfrac{1}{2} \cdot (-2)\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        -1 - \dfrac{1}{2} \cdot 1\\
      \end{array}
    \right) \\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        0 & - \dfrac{5}{2} & 4\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        - \dfrac{3}{2}\\
      \end{array}
    \right)
  \end{eqnarray}
}

となる．次に，2 列目に対して，計算を行う．ここでは， {\dfrac{(3, 2)\mbox{成分}}{(2, 2)\mbox{成分}}} の値を計算し，3 行目 から， {\dfrac{(3, 2)\mbox{成分}}{(2, 2)\mbox{成分}} \times \mbox{2 行目}} を引く．ここでは， {\dfrac{(3, 2)\mbox{成分}}{(2, 2)\mbox{成分}} = -5} となるので，結果は，

 { \displaystyle
  \begin{eqnarray}
    & & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        0 & - \dfrac{5}{2} & 4\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        - \dfrac{3}{2}\\
      \end{array}
    \right)\\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        0 & - \dfrac{5}{2} + 5 \cdot \dfrac{1}{2} & 4 + 5 \cdot 0\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        - \dfrac{3}{2} + 5 \cdot \dfrac{7}{2}\\
      \end{array}
    \right) \\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & \dfrac{1}{2} & 0\\
        0 & 0 & 4\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        \dfrac{7}{2}\\
        16\\
      \end{array}
    \right) \\
    & \Leftrightarrow & \left(
      \begin{array}{ccc}
        2 & 1 & -2\\
        0 & 1 & 0\\
        0 & 0 & 4\\
      \end{array}
    \right)\left(
      \begin{array}{c}
        x_{1}\\
        x_{2}\\
        x_{3}\\
      \end{array}
    \right) = \left(
      \begin{array}{c}
        1\\
        7\\
        16\\
      \end{array}
    \right)
  \end{eqnarray}
}

となる．最後の処理は必要ないが，ここに書く上で，分数は記述が面倒であるし，見づらくなるので，2 行目だけ，両辺を 2 倍しておく．これは，両辺に左から

 { \displaystyle
  \begin{equation}
    \left(
    \begin{array}{ccc}
      1 & 0 & 0\\
      0 & 2 & 0\\
      0 & 0 & 1\\
    \end{array}
    \right)
  \end{equation}
}

をかけることと等しい．これで，前進消去が行えた．後は，後退代入を行っていけばよい．後退代入は，行列をもとの連立方程式の形に直したほうが分かりやすい．連立方程式の形に直すと，

 { \displaystyle
  \begin{equation}
  \left\{
    \begin{array}{cccc}
      2x_{1} & + x_{2} & - 2x_{2} & = 1\\
      & x_{2} & + 0x_{3} & = 7\\
      & & 4x_{3} & = 16\\
    \end{array}
  \right.
  \end{equation}
}

となる．ここで，あえて， {0x_{3}} を入れてある．この理由は，後々分かる．さて，後退代入は， {x_{3}, x_{2}, x_{1}} の順に計算してゆく．方程式を見ると， {x_{3} = 4} となることが直ちに分かる．次に， {x_{2}} は，以下のように式変形して求める．

 { \displaystyle
  \begin{eqnarray}
    & & x_{2} + 0 \cdot x_{3} = 7\\
    & \Leftrightarrow & x_{2} = 7 - 0 \cdot x_{3}\\
    & \Leftrightarrow & x_{2} = 7 - 0 \cdot 4\\
    & \therefore & x_{2} = 7
  \end{eqnarray}
}

この場合は， {x_{3}} の係数が 0 であったため，あまり必要性を感じないかも知れない．なので，先程と同様に {x_{1}} を計算すると，

 { \displaystyle
  \begin{eqnarray}
    & & 2x_{1} + x_{2} - 2x_{3} = 1\\
    & \Leftrightarrow & 2x_{1} = 1 - x_{2} + 2x_{3}\\
    & \Leftrightarrow & 2x_{1} = 1 - 7 + 2 \cdot 4\\
    & \Leftrightarrow & 2x_{1} = 2\\
    & \therefore & x_{1} = 1
  \end{eqnarray}
}

となる．これで，すべての解が計算できた．

アルゴリズム とサンプルプログラム

これのアルゴリズムは以下のようになる．

しかし，これでは，対角成分が 0 の場合は，計算が出来ない．そこで，ピボット選択ということを行う．これは，現在の列のうち，絶対値が最大となるものを探し，行を交換するというものである．これは，対角成分が 0 であるときの対処法だけではなく，丸め誤差の影響が小さくなる．本当は完全ピボット選択の方が良いが，プログラムを組むのが面倒なので，部分ピボット選択を行うガウスの消去法のプログラムを以下に示す．

このプログラムは，行列を 1 次元配列で確保しているので，row 行 col の行列 {A} の {(i, j)} 成分 {a_{ij}} は，A[i * col + j] となる．

参考文献

ガウスの消去法 - Wikipedia