あもんノート TOP へ

数値計算 for VBA

Excel は Microsoft 社製の表計算ソフトで、 多くの Windows パソコンにプレインストールされています。 グラフの作成や図のドローイング、 さらに内蔵する VBA (Visual BASIC for Applications) によりプログラミングもできるので、 日常的業務のみならず数理分野においても重宝します。 簡易なワープロやデータベースとしても使うことができるため、 パソコン世界における「何でも屋」と呼べる定番ソフトです。

ただし VBA には 複素数 の型が用意されていません。 複素数を含んだ数値計算を行うには、それを 構造体 (ユーザー定義型) で定義してやる必要があります。 ここでは構造体により複素数および複素行列の型を作ってやり、 さらに関連した関数を用意し、数値計算に利用することにします。

以下はその基礎となるモジュールです。 全体をコピーして VBA の標準モジュールに貼り付ければ使用できます。 VBA の開発環境は Excel のメニューから「開発→Visual BASIC」で呼び出せます。 標準モジュールはさらに「挿入→標準モジュール」により呼び出せます。

'=== 複素数と複素行列の型 ===

Type complex
        re As Double '実部
        im As Double '虚部
End Type

Type matrix
        row As Long '行数
        col As Long '列数
        ele(9, 9) As complex '要素
End Type

'=== 複素数の関数 ===

'複素数の実表示
Function cn(a As Double, b As Double) As complex
        cn.re = a
        cn.im = b
End Function

'複素数の絶対値
Function cabs(a As complex) As Double
        cabs = (a.re * a.re + a.im * a.im) ^ 0.5
End Function

'複素共役
Function cconj(a As complex) As complex
        cconj.re = a.re
        cconj.im = -a.im
End Function

'複素数の和
Function csum(a As complex, b As complex) As complex
        csum.re = a.re + b.re
        csum.im = a.im + b.im
End Function

'複素数の差
Function cdif(a As complex, b As complex) As complex
        cdif.re = a.re - b.re
        cdif.im = a.im - b.im
End Function

'複素数の積
Function cpro(a As complex, b As complex) As complex
        cpro.re = a.re * b.re - a.im * b.im
        cpro.im = a.re * b.im + a.im * b.re
End Function

'複素数の除算
Function cdiv(a As complex, b As complex) As complex
        cdiv = cpro(cpro(a, cconj(b)), cn(cabs(b) ^ -2, 0))
End Function

'複素数の指数関数
Function cexp(a As complex) As complex
        cexp.re = Exp(a.re) * Cos(a.im)
        cexp.im = Exp(a.re) * Sin(a.im)
End Function

'複素数の文字列化
Function ctos(a As complex) As String
        ctos = "(" + Str(a.re) + "," + Str(a.im) + ")"
End Function

'文字列の複素数化
Function stoc(s As String) As complex
        Dim t As String: t = Trim(s)
        If InStr(t, ",") < 1 Then
                stoc = cn(Val(t), 0)
        Else
                t = Mid(t, 2, Len(t) - 2)
                Dim p As Long: p = InStr(t, ",")
                stoc = cn(Val(Left(t, p - 1)), Val(Mid(t, p + 1)))
        End If
End Function

'=== 複素行列の関数 ===

'零行列
Function zero(row As Long, col As Long) As matrix
        zero.row = row: zero.col = col
        Dim i As Long, j As Long
        For i = 1 To row: For j = 1 To col
                zero.ele(i, j) = cn(0, 0)
        Next: Next
End Function

'単位行列
Function identity(d As Long) As matrix
        Dim i As Long, s As matrix
        s = zero(d, d)
        For i = 1 To d
                s.ele(i, i) = cn(1, 0)
        Next
        identity = s
End Function

'複素数と行列の積
Function cmpro(c As complex, a As matrix) As matrix
        cmpro.row = a.row: cmpro.col = a.col
        Dim i As Long, j As Long
        For i = 1 To a.row: For j = 1 To a.col
                cmpro.ele(i, j) = cpro(c, a.ele(i, j))
        Next: Next
End Function

'行列の和
Function msum(a As matrix, b As matrix) As matrix
        msum.row = a.row: msum.col = a.col
        Dim i As Long, j As Long
        For i = 1 To a.row: For j = 1 To a.col
                msum.ele(i, j) = csum(a.ele(i, j), b.ele(i, j))
        Next: Next
End Function

'行列の差
Function mdif(a As matrix, b As matrix) As matrix
        mdif = msum(a, cmpro(cn(-1, 0), b))
End Function

'行列の積
Function mpro(a As matrix, b As matrix) As matrix
        mpro.row = a.row: mpro.col = b.col
        Dim i As Long, j As Long, k As Long, s As complex
        For i = 1 To a.row: For j = 1 To b.col
                s = cn(0, 0)
                For k = 1 To a.col
                        s = csum(s, cpro(a.ele(i, k), b.ele(k, j)))
                Next
                mpro.ele(i, j) = s
        Next: Next
End Function

'行列の複素共役
Function cconjm(a As matrix) As matrix
        cconjm.row = a.row: cconjm.col = a.col
        Dim i As Long, j As Long
        For i = 1 To a.row: For j = 1 To a.col
                cconjm.ele(i, j) = cconj(a.ele(i, j))
        Next: Next
End Function

'転置行列
Function trans(a As matrix) As matrix
        trans.row = a.col: trans.col = a.row
        Dim i As Long, j As Long
        For i = 1 To a.row: For j = 1 To a.col
                trans.ele(j, i) = a.ele(i, j)
        Next: Next
End Function

'エルミート共役
Function hermit(a As matrix) As matrix
        hermit = trans(cconjm(a))
End Function

'トレース
Function trace(a As matrix) As complex
        Dim i As Long, s As complex
        s = cn(0, 0)
        For i = 1 To a.row
                s = csum(s, a.ele(i, i))
        Next
        trace = s
End Function

'交換子
Function commutator(a As matrix, b As matrix) As matrix
        commutator = mdif(mpro(a, b), mpro(b, a))
End Function

'反交換子
Function anticommutator(a As matrix, b As matrix) As matrix
        anticommutator = msum(mpro(a, b), mpro(b, a))
End Function

'=== 行列式と余因子 ===

'行列式
Function det(a As matrix) As complex
        Select Case a.row
        Case 1
                det = a.ele(1, 1)
        Case Else
                Dim j As Long, s As complex
                s = cn(0, 0)
                For j = 1 To a.row
                        s = csum(s, cpro(a.ele(1, j), cofactor(a, 1, j)))
                Next
                det = s
        End Select
End Function

'余因子
Function cofactor(a As matrix, i As Long, j As Long) As complex
        Dim b As matrix
        b.row = a.row - 1: b.col = b.row
        Dim x As Long, y As Long, xd As Long, yd As Long
        For x = 1 To b.row: For y = 1 To b.col
                xd = x: If x >= i Then xd = x + 1
                yd = y: If y >= j Then yd = y + 1
                b.ele(x, y) = a.ele(xd, yd)
        Next: Next
        cofactor = cpro(cn((-1) ^ (i + j), 0), det(b))
End Function

'余因子行列
Function comatrix(a As matrix) As matrix
        comatrix.row = a.row: comatrix.col = a.col
        Dim i As Long, j As Long
        For i = 1 To a.row: For j = 1 To a.row
                comatrix.ele(i, j) = cofactor(a, j, i)
        Next: Next
End Function
大文字で始まる単語はシステムの予約語です。 Type は構造体、Function は関数(返り値のあるメソッド)、 Sub はサブルーチン(返り値のないメソッド)の宣言を意味します。 Dim は変数の宣言、Long は長整数型、Double は倍精度浮動小数点数型、 String は文字列型を意味します。 上のモジュールにより、さらに complex が複素数型、 matrix が複素行列型を意味するようになり、 これらの型に関連した様々な関数を利用できるようになるというわけです。

以下は行列式と余因子行列の計算プログラムです。 ワークシートのコードウィンドウに貼り付ければ使用できます。 使い方ですが、対応するワークシート左上から正方行列の各成分を書き並べます。 ただし複素数 x+iy は (x,y) と記します。 実数なら普通の数値の記法で大丈夫です。 そうして det_comatrix_out を実行すると行列式と余因子行列が同じワークシード上に出力されます。 実行は Excel のメニューの「開発→マクロ」から行います。 (一般にマクロ実行のイベントは Excel 上の図形やコントロールに貼り付け可能です。)

'行列式と余因子行列の出力
Sub det_comatrix_out()
        Dim i As Long, j As Long, n As Long
        Dim a As matrix, b As matrix
        n = 0
        Do
                n = n + 1
        Loop Until Cells(n, 1).Value = ""
        n = n - 1
        If n > 1 Then
                a.row = n: a.col = n
                For i = 1 To n: For j = 1 To n
                        a.ele(i, j) = stoc(Cells(i, j).Value)
                Next: Next
                Cells(n + 2, 1).Value = ctos(det(a))
                b = comatrix(a)
                For i = 1 To n: For j = 1 To n
                        Cells(n + 3 + i, j).Value = ctos(b.ele(i, j))
                Next: Next
        End If
End Sub
以下は SU(3) の構造定数の計算プログラムです。 ワークシートのコードウィンドウに貼り付ければ使用できます。 (連続群論入門の章参照。)

'ゲルマン行列
Function gellmann(a As Long) As matrix
        Dim g As matrix: g = zero(3, 3)
        Select Case a
        Case 1
                g.ele(1, 2) = cn(1, 0): g.ele(2, 1) = cn(1, 0)
        Case 2
                g.ele(1, 2) = cn(0, -1): g.ele(2, 1) = cn(0, 1)
        Case 3
                g.ele(1, 1) = cn(1, 0): g.ele(2, 2) = cn(-1, 0)
        Case 4
                g.ele(1, 3) = cn(1, 0): g.ele(3, 1) = cn(1, 0)
        Case 5
                g.ele(1, 3) = cn(0, -1): g.ele(3, 1) = cn(0, 1)
        Case 6
                g.ele(2, 3) = cn(1, 0): g.ele(3, 2) = cn(1, 0)
        Case 7
                g.ele(2, 3) = cn(0, -1): g.ele(3, 2) = cn(0, 1)
        Case 8
                g.ele(1, 1) = cn(1, 0): g.ele(2, 2) = cn(1, 0): g.ele(3, 3) = cn(-2, 0)
                g = cmpro(cn(3 ^ (-0.5), 0), g)
        End Select
        gellmann = g
End Function

'SU(3)構造定数出力
Sub su3structure()
        Dim i As Long: i = 1
        Dim a As Long, b As Long, c As Long
        Dim f As complex, m As matrix
        For a = 1 To 8: For b = a + 1 To 8: For c = b + 1 To 8
                m = mpro(commutator(gellmann(a), gellmann(b)), gellmann(c))
                f = cpro(cn(0, -1 / 4), trace(m))
                If Not (f.re = 0 And f.im = 0) Then
                        i = i + 1
                        Cells(i, 1).Value = Str(a) + Str(b) + Str(c)
                        Cells(i, 2).Value = ctos(f)
                End If
        Next: Next: Next
End Sub
以下は電子陽電子対消滅の断面積数値計算のコードです。 ワークシートのコードウィンドウに貼り付ければ使用できます。 (素粒子論の計算の章参照。)

'円周率
Const pi = 3.14159265358979

'ガンマ行列(下付き添字)
Function gamma(a As Long) As matrix
        Dim g As matrix: g = zero(4, 4)
        Select Case a
        Case 0
                g.ele(1, 3) = cn(1, 0): g.ele(2, 4) = cn(1, 0)
                g.ele(3, 1) = cn(1, 0): g.ele(4, 2) = cn(1, 0)
        Case 1
                g.ele(1, 4) = cn(1, 0): g.ele(2, 3) = cn(1, 0)
                g.ele(3, 2) = cn(-1, 0): g.ele(4, 1) = cn(-1, 0)
        Case 2
                g.ele(1, 4) = cn(0, -1): g.ele(2, 3) = cn(0, 1)
                g.ele(3, 2) = cn(0, 1): g.ele(4, 1) = cn(0, -1)
        Case 3
                g.ele(1, 3) = cn(1, 0): g.ele(2, 4) = cn(-1, 0)
                g.ele(3, 1) = cn(-1, 0): g.ele(4, 2) = cn(1, 0)
        Case 5
                g.ele(1, 1) = cn(1, 0): g.ele(2, 2) = cn(1, 0)
                g.ele(3, 3) = cn(-1, 0): g.ele(4, 4) = cn(-1, 0)
        End Select
        gamma = g
End Function

'ディラック粒子の波動関数 ( s はヘリシティ*2)
Function dirac_u(s As Long, p As Double, th As Double) As matrix
        Dim e As Double: e = (p ^ 2 + 1) ^ 0.5
        Dim rp As Double: rp = (e + p) ^ 0.5
        Dim rm As Double: rm = (e - p) ^ 0.5
        Dim sth As Double: sth = Sin(th / 2)
        Dim cth As Double: cth = Cos(th / 2)
        Dim a As matrix: a = zero(4, 1)
        Select Case s
        Case 1
                a.ele(1, 1) = cn(rp * cth, 0): a.ele(2, 1) = cn(rp * sth, 0)
                a.ele(3, 1) = cn(rm * cth, 0): a.ele(4, 1) = cn(rm * sth, 0)
        Case -1
                a.ele(1, 1) = cn(0, -rm * sth): a.ele(2, 1) = cn(0, rm * cth)
                a.ele(3, 1) = cn(0, -rp * sth): a.ele(4, 1) = cn(0, rp * cth)
        End Select
        dirac_u = a
End Function

'ディラック反粒子の波動関数
Function dirac_v(s As Long, p As Double, th As Double) As matrix
        dirac_v = cmpro(cn(0, -1), mpro(gamma(2), cconjm(dirac_u(s, p, th))))
End Function

'ディラック共役
Function dirac_bar(a As matrix) As matrix
        dirac_bar = mpro(hermit(a), gamma(0))
End Function

'光子の波動関数 ( m は時空添字、l はヘリシティ)
Function photon_eps(m As Long, l As Long, th As Double) As complex
        Select Case m
        Case 0: photon_eps = cn(0, 0)
        Case 1: photon_eps = cn(2 ^ (-0.5) * Cos(th), 0)
        Case 2: photon_eps = cn(0, 2 ^ (-0.5) * l)
        Case 3: photon_eps = cn(2 ^ (-0.5) * (-Sin(th)), 0)
        End Select
End Function

'光子の波動関数*ガンマ行列 ( l はヘリシティ)
Function eps_sla(l As Long, th As Double) As matrix
        Dim a As matrix, m As Long
        a = zero(4, 4)
        For m = 1 To 3
                a = msum(a, cmpro(photon_eps(m, l, th), gamma(m)))
        Next
        eps_sla = a
End Function

'直接グラフ A/α
Function graph_a( _
s As Long, sd As Long, l As Long, ld As Long, p As Double, th As Double _
) As complex
        Dim e As Double: e = (p ^ 2 + 1) ^ 0.5
        Dim a As matrix: a = identity(4)
        a = msum(a, cmpro(cn(-e * Sin(th), 0), gamma(1)))
        a = msum(a, cmpro(cn(p - e * Cos(th), 0), gamma(3)))
        a = mpro(mpro(eps_sla(-ld, th + pi), a), eps_sla(-l, th))
        a = mpro(mpro(dirac_bar(dirac_v(sd, p, pi)), a), dirac_u(s, p, 0))
        graph_a = cpro(cn(0, 2 * pi / e / (e - p * Cos(th))), trace(a))
End Function

'散乱振幅 M/α
Function s_amp( _
s As Long, sd As Long, l As Long, ld As Long, p As Double, th As Double _
) As complex
        Dim a As complex: a = graph_a(s, sd, l, ld, p, th)
        Dim b As complex: b = graph_a(s, sd, ld, l, p, th + pi)
        s_amp = csum(a, b)
End Function

'微分断面積 dσ/dΩ*(m/α)^2
Function dif_cs( _
s As Long, sd As Long, l As Long, ld As Long, p As Double, th As Double _
) As Double
        Dim e As Double: e = (p ^ 2 + 1) ^ 0.5
        dif_cs = cabs(s_amp(s, sd, l, ld, p, th)) ^ 2 / (256 * pi ^ 2 * p * e)
End Function

'断面積 σ*(m/α)^2
Function cs( _
s As Long, sd As Long, l As Long, ld As Long, p As Double, n As Long _
) As Double
        Dim dth As Double: dth = pi / 2 / n: If l <> ld Then dth = pi / n
        Dim ov As Double, nv As Double, a As Double, i As Long
        a = 0: ov = 0
        For i = 1 To n
                nv = Sin(dth * i) * dif_cs(s, sd, l, ld, p, dth * i)
                a = a + (ov + nv) / 2 * dth
                ov = nv
        Next
        cs = a * 2 * pi
End Function

'断面積の数値出力
Sub cs_output()
        Dim n As Long: n = 100 'θ積分分割数
        Cells(1, 1).Value = "p/m"
        Cells(1, 2).Value = "++++"
        Cells(1, 3).Value = "+++-"
        Cells(1, 4).Value = "++--"
        Cells(1, 5).Value = "+-++"
        Cells(1, 6).Value = "+-+-"
        Dim i As Long, p As Double
        For i = 0 To 60
                p = i / 10: If i = 0 Then p = 0.001
                Cells(i + 2, 1).Value = p
                Cells(i + 2, 2).Value = cs(1, 1, 1, 1, p, n)
                Cells(i + 2, 3).Value = cs(1, 1, 1, -1, p, n)
                Cells(i + 2, 4).Value = cs(1, 1, -1, -1, p, n)
                Cells(i + 2, 5).Value = cs(1, -1, 1, 1, p, n)
                Cells(i + 2, 6).Value = cs(1, -1, 1, -1, p, n)
        Next
End Sub
以下はゼータ関数数値出力のためのコードです。 ワークシートのコードウィンドウに貼り付ければ使用できます。 (ゼータ関数の章参照。)

'円周率
Const pi = 3.14159265358979

'実数の複素べき
Function cpow(t As Double, z As complex) As complex
        cpow = cexp(cpro(z, cn(Log(t), 0)))
End Function

'複素数の三角関数
Function csin(z As complex) As complex
        Dim a As complex
        a = cdif(cexp(cpro(cn(0, 1), z)), cexp(cpro(cn(0, -1), z)))
        csin = cpro(cn(0, -0.5), a)
End Function

'被積分関数
Function itgr1(z As complex, t As Double) As complex
        itgr1 = cpro(cpow(t, z), cn(1 / t / Exp(t), 0))
End Function

Function itgr2(z As complex, th As Double) As complex
        Dim a As complex: a = cexp(cpro(z, cn(0, th - pi)))
        Dim b As complex: b = cdif(cexp(cexp(cn(0, th))), cn(1, 0))
        itgr2 = cdiv(a, b)
End Function

Function itgr3(z As complex, t As Double) As complex
        itgr3 = cpro(cpow(t, z), cn(1 / t / (Exp(t) - 1), 0))
End Function

'ガンマ関数 ( tmax は積分上限値、n は積分分割数)
Function gamma_f(z As complex, tmax As Double, n As Long) As complex
        If z.re < 0.5 Then
                Dim a As complex: a = cdiv(cn(pi, 0), csin(cpro(cn(pi, 0), z)))
                gamma_f = cdiv(a, gamma_f(cdif(cn(1, 0), z), tmax, n))
        Else
                Dim dt As Double: dt = tmax / n
                Dim i As Long, t1 As Double, t2 As Double
                Dim s As complex: s = cdiv(cpow(dt, z), z)
                For i = 1 To n - 1
                        t1 = dt * i: t2 = dt * (i + 1)
                        s = csum(s, cpro(csum(itgr1(z, t1), itgr1(z, t2)), cn(dt / 2, 0)))
                Next
                gamma_f = s
        End If
End Function

'ゼータ関数 ( tmax は積分上限値、n は積分分割数)
Function zeta_f(z As complex, tmax As Double, n As Long) As complex
        If z.re > 0 And z.im = 0 And Int(z.re) = z.re Then
                Dim eps As Double: eps = 0.0001
                Dim a As complex: a = zeta_f(csum(z, cn(0, eps)), tmax, n)
                Dim b As complex: b = zeta_f(cdif(z, cn(0, eps)), tmax, n)
                zeta_f = cpro(csum(a, b), cn(0.5, 0))
        Else
                Dim i As Long, t1 As Double, t2 As Double
                Dim dth As Double: dth = 2 * pi / n
                Dim s1 As complex: s1 = cn(0, 0)
                For i = 0 To n - 1
                        t1 = i * dth: t2 = (i + 1) * dth
                        s1 = csum(s1, cpro(csum(itgr2(z, t1), itgr2(z, t2)), cn(dth / 2, 0)))
                Next
                s1 = cpro(s1, cn(1 / 2 / pi, 0))
                Dim dt As Double: dt = (tmax - 1) / n
                Dim s2 As complex: s2 = cn(0, 0)
                For i = 0 To n - 1
                        t1 = 1 + i * dt: t2 = 1 + (i + 1) * dt
                        s2 = csum(s2, cpro(csum(itgr3(z, t1), itgr3(z, t2)), cn(dt / 2, 0)))
                Next
                s2 = cpro(s2, cpro(csin(cpro(cn(pi, 0), z)), cn(1 / pi, 0)))
                zeta_f = cpro(gamma_f(cdif(cn(1, 0), z), tmax, n), csum(s1, s2))
        End If
End Function

'実軸上のゼータ関数数値出力
Sub zeta_f_output()
        Dim i As Long, x As Double, y As complex
        For i = 0 To 100
                x = -5 + i / 10
                If Not x = 1 Then
                        y = zeta_f(cn(x, 0), 50, 10000)
                        Cells(i + 1, 1).Value = x
                        Cells(i + 1, 2).Value = y.re
                End If
        Next
End Sub
プログラミング言語にある程度精通しているが VBA は知らないという方へ。

VBA は基本的には静的型付けの言語です。 ただし型宣言をしないと Variant 型と呼ばれる全てのプリミティブを内包する型になり、 その扱いはインタプリタにゆだねられます。 すなわち動的型付けっぽい一面も有しています。 ただし VBA の性質を考えると、型宣言はちゃんとした方が良いでしょう。

関数にいわゆる return 文がないのは特徴的でしょう。 関数と同名の変数に値をセットして関数のブロック(プロシージャ)を終えると、 その値が返されます。 大文字と小文字は区別されません。 予約語の先頭文字は開発環境が自動的に大文字にします。 Cells は Excel のセルオブジェクトで、 この他 Excel や OS が提供する様々なオブジェクトを VBA から操作できます。 クラスの機能も用意されていますが、 メンバ関数の名前に予約語が使えないこともあり、 これは使い物にならないと考えた方が良いでしょう。 一方でご覧のように構造体が用意されていて、 数値計算や事務処理には十分な機能を有しています。

コードを見てもわかるように、 BASIC 系の言語の特徴は構文が英語に近いことで、 C 系の言語のように括弧が多用されることはありません。 このため初心者が馴染みやすい言語と考えられています。

あもんノート TOP へ
inserted by FC2 system