あもんノート 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): Exit Function
    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 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
    For i = 1 To 9
        If Cells(i, 1).Value = "" Then Exit For
    Next
    n = i - 1: If n < 2 Then Exit Sub
    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 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))
        Exit Function
    End If
    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 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))
        Exit Function
    End If
    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 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 型と呼ばれる全てのプリミティブを内包する型になり、 その扱いはインタプリタにゆだねられます。 すなわち動的型付けっぽい一面も有しています。 型宣言を強要する場合はコードの最初に Option Explicit と書きます。

関数にいわゆる return 文がないのは特徴的でしょう。 関数と同名の変数に値をセットして関数のブロック(プロシージャ)を終えると、 その値が返されます。 強制的にブロックを終えるコマンド(キーワード)に Exit があります。

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

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