あもんノート 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 へ