連立方程式(1) ガウスの消去法

連立方程式を解く方法で、大変広く知られているアルゴリズムです。何も知らずに連立方程式をPCに解かせようと思うと、プログラムを考えたり、ループを考えたりして何時間、何日もかかってしまいますが、一度このアルゴリズムを覚えてしまうと数分で書けますので、是非覚えておきたいコードです。
ここで紹介するプログラムは本当にシンプルに書いてありますが、欠点は係数a中で対角行列に0が含まれると計算できなくなることです(0の除算が含まれてしまう)。これは、あらかじめ行列の順序を入れ替えることで回避できます(解法その2参照)。

エクセルでは図のように、係数となるa行列、及び解となるb行列を用意し、この値を読むところから始まります。 サンプルで用意したプログラムでは5変数まで求められるようになっていますが、行列を大きくすれば、いくらでも大きな連立方程式を解くことができます。

ロジック
ここでは、
\( ax=b \)
を解きます。変数の数をnとし成分を明示すると、

\(
\begin{pmatrix}
   a_{11} & \ldots & a_{1n} \\
   \vdots & \ddots & \vdots \\
   a_{n1} & \ldots & a_{nn}
\end{pmatrix}.
\begin{pmatrix}
   x_{1} \\   \vdots \\   x_{n}
\end{pmatrix}
=
\begin{pmatrix}
   b_{1} \\   \vdots \\   b_{n}
\end{pmatrix}
\)
となります。この行列を、まず逆▲行列にします。即ち、次のように対角より下の成分をすべて0にします。プログラム中では§2にあたります。

\(
\begin{pmatrix}
   a_{11} & a_{12}   & \ldots  & a_{1n} \\
   0       & a_{22}   & \ldots & a_{2n} \\
   0       & 0         & \ddots & \vdots \\
   0       &  \ldots  &    0      & a_{nn}
\end{pmatrix}.
\begin{pmatrix}
   x_{1} \\   x_{2} \\ \vdots \\   x_{n}
\end{pmatrix}
=
\begin{pmatrix}
   b_{1} \\  b_{2} \\ \vdots \\   b_{n}
\end{pmatrix}
\)
となります。
行列がこの状態になると、まず\( x_{n}\)が、

\(
x_{n} = b_n/a_{nn}
\)
として求められます。その\( x_{n}\)を使って、\( x_{n-1}\)から\( x_{1}\)に向かって順に解を求めていきます。(§3)





下にプログラムを示します。

Public Const m = 10

Sub test()
  Dim i%, j%, k%, n%
  Dim a(m, m) As Double, b(m) As Double, x(m) As Double

  n = Cells(2, 3)  ' 変数の数を読み込む

  '§1. データの読み込み ------------------------
  For i = 1 To n
    For j = 1 To n
      a(i, j) = Cells(i + 5, j + 1)
    Next j
    b(i) = Cells(i + 5, 8)
  Next i

  '§2. 上三角(▼)行列にする --------------------
  For k = 1 To n
    For i = k + 1 To n
        AI = a(i, k) / a(k, k)
        For j = 1 To n
          a(i, j) = a(i, j) - AI * a(k, j)
        Next j
        b(i) = b(i) - AI * b(k)
    Next i
  Next k

  '§3. 成分nから成分1に向かって解を求めていく ---
  x(n) = b(n) / a(n, n)
  For k = n - 1 To 1 Step -1
    For j = k + 1 To n
      b(k) = b(k) - a(k, j) * x(j)
    Next j
    x(k) = b(k) / a(k, k)
  Next k

  '§4. シートへ結果を書き出す ------------------
  For i = 1 To n
    For j = 1 To n
      Cells(i + 12, j + 1) = a(i, j)
    Next j
    Cells(i + 12, 8) = b(i)
    Cells(i + 12, 9) = x(i)
  Next i

End Sub








ダウンロード

0 件のコメント:

コメントを投稿