行列計算に関するExcelマクロ

設計シートで用いる行列計算に関するマクロのプログラム(関数)を紹介する

関数の説明
サンプルコード
まとめ

関数の説明

アルゴリズムに関する具体的な解説は他にもっと詳しいサイトがあると思うのでそちらに譲り，ここではどんな引数や返り値を使用するかの説明にとどめる

連立方程式のソルバー

連立一次方程式\({\bf Ax=b}\)のソルバーはガウスの消去法を用いる

直接法だとか間接法だとか，計算速度がNの2乗に比例するだとかいろいろあるがここでは説明しない

Public Function ans(ByVal N As Integer, ByRef A() As Double, ByRef b() As Double) As Double()

関数名はans．引数は行列の要素数n，n×n行列A，n×1行列b．戻り値はn×1行列x

行列の要素数は，例えばA(0:n-1,0:n-1)で宣言した場合にnであることに注意する

逆行列の計算

逆行列の計算には掃き出し法を用いる

Public Function inverse(ByVal N As Integer, ByRef b() As Double) As Double()

関数名はinverse，引数は行列の要素数n，n×n行列b．戻り値はbの逆行列

行列の掛け算

行列の掛け算は

Public Function mult(m1() As Double, m2() As Double) As Double()

関数名はmult，引数は掛け算したい行列m1とm2．戻り値はm1とm2の積．

このコードは内部で行列の要素数を計算するようになっているので行列の要素数を引数にする必要はない

転置行列の計算

転置行列の計算は

Public Function Transpose(ByVal M As Integer, ByVal N As Integer, ByRef A() As Double) As Double()

関数名はTranspose，引数は行列の要素数(m,n)とm×n行列A．戻り値はAの転置行列

サンプルコード

自分で書いたものやネットで拾ってきたものも混ざっているので使用は自己責任で

ガウスの消去法

Public Function ans(ByVal N As Integer, ByRef A() As Double, ByRef b() As Double) As Double()
'ガウスの消去法によって，連立方程式Ax=bを解く．
'Aはn×n行列，x，bはn×1行列
'参考文献：「解析塾秘伝！有限要素法のつくり方」　p.82-

'変数の宣言
Dim pivot As Double 'マトリックスの体格成分
Dim p As Double '計算に使用するマトリックスの成分
Dim x() As Double '解を入れる配列
ReDim x(N - 1)

'カウンター
Dim r As Integer
Dim c As Integer
Dim rr As Integer
Dim cc As Integer

'画面更新の非表示
Application.ScreenUpdating = False


'前進消去
For r = 0 To N - 1
    pivot = A(r, r) '対角成分をpivotに代入
    For c = r To N - 1
        A(r, c) = A(r, c) / pivot
    Next c
    b(r) = b(r) / pivot
    For rr = r + 1 To N - 1
        p = A(rr, r)
        For cc = r To N - 1
            A(rr, cc) = A(rr, cc) - p * A(r, cc)
        Next cc
        b(rr) = b(rr) - p * b(r)
    Next rr
Next r

'後退代入
For r = N - 1 To 0 Step -1
    x(r) = b(r)
    For c = r + 1 To N - 1
        x(r) = x(r) - A(r, c) * x(c)
    Next c
Next r

ans = x

End Function

逆行列の計算

Public Function inverse(ByVal N As Integer, ByRef b() As Double) As Double()
'掃き出し法により逆行列を計算する
'参考資料：https://www.econ.nagoya-cu.ac.jp/~kamiyama/siryou/minv.html

Dim nn As Integer
Dim p As Double 'ピボット
Dim c As Double '行
Dim x As Double
Dim y As Double
Dim A() As Double
Dim Binv() As Double

nn = N + N

ReDim A(N - 1, nn - 1)
ReDim Binv(N - 1, N - 1)

'左半分に逆行列を求めたい行列を入れる
For i = 0 To N - 1
    For j = 0 To N - 1
        A(i, j) = b(i, j)
    Next j
Next i

'行列の右半分に単位行列を入れる
For i = 0 To N - 1
    For j = N To nn - 1
        A(i, j) = 0
    Next j
    A(i, i + N) = 1
Next i

'掃き出し法の基本操作
For k = 0 To N - 1
    p = 0: c = 0
    'k行以降で絶対値が最も大きい行を調べる
    For i = k To N - 1
        If Abs(A(i, k)) > p Then
            p = Abs(A(i, k))
            c = i
        End If
    Next i
    'k行目とc行目の成分を入れ替える
    If c <> k Then
        For j = k To nn - 1
            y = A(k, j)
            A(k, j) = A(c, j)
            A(c, j) = y
        Next j
    End If
    'k行目以降の対角成分をk行目の対角成分で割る
    x = A(k, k)
    For j = k To nn - 1
        A(k, j) = A(k, j) / x
    Next j
    For i = 0 To N - 1
        If i <> k Then
            x = A(i, k)
            For j = k To nn - 1
                A(i, j) = A(i, j) - x * A(k, j)
            Next j
        End If
    Next i
Next k

For i = 0 To N - 1
    For j = N To nn - 1
        Binv(i, j - N) = A(i, j)
    Next j
Next i

inverse = Binv

End Function

行列の掛け算

Public Function mult(m1() As Double, m2() As Double) As Double()
'm1[l,m]*m2[m,n]の行列の掛け算を行う

Dim L As Integer
Dim M As Integer
Dim N As Integer
Dim sum As Double
Dim ans() As Double

'カウンター
Dim i As Integer
Dim j As Integer
Dim k As Integer

'行列の大きさを取得
L = UBound(m1, 1)
M = UBound(m1, 2)
N = UBound(m2, 2)

'配列の大きさを定義
ReDim ans(L, N)

For i = 0 To L
    For j = 0 To N
    sum = 0
     For k = 0 To M
        sum = sum + m1(i, k) * m2(k, j)
    Next k
    ans(i, j) = sum
    Next j
Next i

mult = ans

End Function

転置行列の計算

Public Function Transpose(ByVal M As Integer, ByVal N As Integer, ByRef A() As Double) As Double()
'転置行列を計算する
'Aはm×n行列，Atはn×m行列

'変数の宣言
Dim At() As Double '解を入れる配列
ReDim At(N - 1, M - 1) '解を入れる配列

'カウンター
Dim r As Integer
Dim c As Integer

For r = 0 To M - 1
    For c = 0 To N - 1
        At(c, r) = A(r, c)
    Next c
Next r

Transpose = At

End Function