Excel VBA: 連立方程式

連立方程式を解く方法その2（その1はこちら）で、やはりこれも大変広く知られているアルゴリズムです。その1で紹介したガウスの消去法では対角行列に0成分が含まれると計算できなくなりました（もちろん少し追加すれば良いのですが）。このプログラムではそれを回避するアルゴリズムを入れてあります。
行列aの中から絶対値が最大となる成分を探し、それを対角に持ってくるようなアルゴリズムです。これを入れることで、たとえば\(x_2=x_3\)のように解が一意にならないような場合を除き、どんな連立方程式でも解くことができます。

さて、エクセルでは図のように係数となる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}
\)
となります。

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

Public Const m = 10

Sub test()
Dim i%, j%, k%, n%, np1%
Dim irow%, jcol%
Dim a(m, m) As Double
Dim amax#, tmp#
Dim ID(m) As Integer

'1. データの読み込み
'-----------------------------------------------------
'項数を取得
n = Cells(2, 3)
np1 = n + 1 '定数項bをaのn+1列に格納する

For i = 1 To n
For j = 1 To n
a(i, j) = Cells(i + 5, j + 1)
Next j
a(i, np1) = Cells(i + 5, 8)
ID(i) = 0 '計算終了行フラグ
Next i

'2. 計算
'-----------------------------------------------------
For k = 1 To n
'2-1 最大値を見つける
amax = 0
For i = 1 To n
If ID(i) = 0 Then
For j = 1 To n
If Abs(a(i, j)) > amax Then
amax = a(i, j) '最大値を格納
irow = i 'その時の行
jcol = j 'その時の列
End If
Next j
End If
Next i

'2-2 行列の入れ替え
'最大値を対角成分に持ってくる
For j = 1 To np1
tmp = a(irow, j)
a(irow, j) = a(jcol, j)
a(jcol, j) = tmp
Next j
ID(jcol) = 1 'その行を計算したよ、フラグ

'2-3 対角成分を1にする
tmp = a(jcol, jcol)
For j = 1 To np1
a(jcol, j) = a(jcol, j) / tmp
Next j

'2-4 対角成分以外を0にする
For i = 1 To n
If i <> jcol Then
tmp = a(i, jcol)
For j = 1 To np1
a(i, j) = a(i, j) - tmp * a(jcol, j)
Next j
End If
Next i

Next k

'3. 解の出力
'-----------------------------------------------------
For i = 1 To n
For j = 1 To n
Cells(i + 12, j + 1) = a(i, j)
Next j
Cells(i + 12, 8) = a(i, np1)
Next i

End Sub

ダウンロード

連立方程式を解く方法で、大変広く知られているアルゴリズムです。何も知らずに連立方程式を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

ダウンロード

Excel VBA

連立方程式(2) ガウス・ジョルダン法

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