-自作４倍精度クラスでの処理時間- 34桁関数モジュール(clstaby_preci34)
ほぼ同一のプログラムで d＝0.0001 t＝1 の条件で処理時間を比較すると
約2125倍遅くなる。
- 結果 -
倍精度時 1分06秒
自作４倍精度クラス使用 35時間57分54秒
（ムーアの法則だと１０年で1000倍なので11年分処理が遅くなる。 注意：改善余地あり。）
◆◆◆◆ ピタゴラス３体問題 ◆◆◆◆ ※別途クラスモジュール：34桁関数モジュールが必要
Option Explicit
Option Base 1
Dim m(3) As String
Public Const G As String = "1.00000000000000000000000000000000000E+0000" '万有引力定数
Sub メインルーチン()
Dim cls34 As clstaby_preci34
Set cls34 = New clstaby_preci34 ' インスタンスを作成
Dim t As String, dt As String ' 時間
Dim rr(3, 3) As String ' rr（識別番号,軸成分）
Dim vv(3, 3) As String ' vv（識別番号,軸成分）
Dim ans As Variant
Dim cnt0 As Long, cnt1 As Long
Dim sh As Worksheet
Dim 終了時刻 As String
Set sh = ThisWorkbook.Worksheets(1)
sh.Range("C10:AD65536").ClearContents ' データ消去
sh.Range("N2") = Time ' 開始時刻
'初期条件
m(1) = "3.00000000000000000000000000000000000E+0000" ' 質量１
m(2) = "4.00000000000000000000000000000000000E+0000" ' 質量２
m(3) = "5.00000000000000000000000000000000000E+0000" ' 質量３
rr(1, 1) = "1.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(1, 2) = "3.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(1, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(2, 1) = "-2.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(2, 2) = "-1.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(2, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(3, 1) = "1.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(3, 2) = "-1.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
rr(3, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(1, 1) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(1, 2) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(1, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(2, 1) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(2, 2) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(2, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(3, 1) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(3, 2) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
vv(3, 3) = "0.00000000000000000000000000000000000E+0000" '（識別番号,軸成分）
終了時刻 = "1.00000000000000000000000000000000000E+0000"
t = "0.00000000000000000000000000000000000E+0000"
dt = "0.00010000000000000000000000000000000E+0000" '0.0001
cnt0 = 0
cnt1 = 0
With sh '表記
.Cells(9, 4).Value = "t" 't
.Cells(9, 5).Value = "x1" 'x1
.Cells(9, 6).Value = "y1" 'y1
.Cells(9, 7).Value = "x2" 'x2
.Cells(9, 8).Value = "y2" 'y2
.Cells(9, 9).Value = "x3" 'x3
.Cells(9, 10).Value = "y3" 'y3
End With
Do While cls34.■34↓15(t) < cls34.■34↓15(終了時刻)
With sh
.Cells(cnt1 + 10, 4).Value = cls34.■34↓15(t) 't
.Cells(cnt1 + 10, 5).Value = cls34.■34↓15(rr(1, 1)) 'Ｍ1 x
.Cells(cnt1 + 10, 6).Value = cls34.■34↓15(rr(1, 2)) 'Ｍ1 y
.Cells(cnt1 + 10, 7).Value = cls34.■34↓15(rr(2, 1)) 'Ｍ2 x
.Cells(cnt1 + 10, 8).Value = cls34.■34↓15(rr(2, 2)) 'Ｍ2 y
.Cells(cnt1 + 10, 9).Value = cls34.■34↓15(rr(3, 1)) 'Ｍ3 x
.Cells(cnt1 + 10, 10).Value = cls34.■34↓15(rr(3, 2)) 'Ｍ3 y
cnt1 = cnt1 + 1
End With
ans = RK4(t, dt, m, rr, vv) ' ルンゲ＝クッタ
rr(1, 1) = ans(1)
rr(1, 2) = ans(2)
rr(1, 3) = ans(3)
vv(1, 1) = ans(4)
vv(1, 2) = ans(5)
vv(1, 3) = ans(6)
rr(2, 1) = ans(7)
rr(2, 2) = ans(8)
rr(2, 3) = ans(9)
vv(2, 1) = ans(10)
vv(2, 2) = ans(11)
vv(2, 3) = ans(12)
rr(3, 1) = ans(13)
rr(3, 2) = ans(14)
rr(3, 3) = ans(15)
vv(3, 1) = ans(16)
vv(3, 2) = ans(17)
vv(3, 3) = ans(18)
cnt0 = cnt0 + 1
t = cls34.■加34(t, dt) ' t = t + dt
Loop
sh.Range("N3") = Time
Set cls34 = Nothing
End Sub
Function RK4(ByVal t, ByVal dt, m, ByVal r0, ByVal v0) As Variant
Dim cls34 As clstaby_preci34
Set cls34 = New clstaby_preci34 ' インスタンスを作成
'下準備--------------------------------------------------------------------------------------
Dim ret(18)
Dim r1(3, 3), r2(3, 3), r3(3, 3), r4(3, 3) ' r(識別番号,軸成分）
Dim v1(3, 3), v2(3, 3), v3(3, 3), v4(3, 3) ' v(識別番号,軸成分）
Dim k1rr(3, 3), k2rr(3, 3), k3rr(3, 3), k4rr(3, 3) ' k1rr(識別番号,軸成分）
Dim k1vv(3, 3), k2vv(3, 3), k3vv(3, 3), k4vv(3, 3) ' k1vv(識別番号,軸成分）
Dim 番号 As Integer, 軸 As Integer
Dim 数字2 As String
Dim buf01 As String
数字2 = "2.00000000000000000000000000000000000E+0000"
For 番号 = 1 To 3
For 軸 = 1 To 3
v1(番号, 軸) = v0(番号, 軸) 'v(識別番号,軸成分) (3,3) 3体×3次元 前処理
r1(番号, 軸) = r0(番号, 軸) 'r(識別番号,軸成分) (3,3) 3体×3次元 前処理
Next 軸
Next 番号
'下準備--------------------------------------------------------------------------------------
'k1
For 番号 = 1 To 3 '識別番号
For 軸 = 1 To 3 '軸成分
k1vv(番号, 軸) = cls34.■乗34(dt, 速度変位成分_func(番号, 軸, t, r1, v1)) '速度変位_総和 m/s^2
k1rr(番号, 軸) = cls34.■乗34(dt, 位置変位成分_func(番号, 軸, t, r1, v1)) '位置変位成分 m/s
Next 軸
Next 番号
'k2
For 番号 = 1 To 3
For 軸 = 1 To 3
v2(番号, 軸) = cls34.■加34(v1(番号, 軸), cls34.■除34(k1vv(番号, 軸), 数字2)) 'v(識別番号,軸成分) (3,3) 前処理
r2(番号, 軸) = cls34.■加34(r1(番号, 軸), cls34.■除34(k1rr(番号, 軸), 数字2)) 'r(識別番号,軸成分) (3,3) 前処理
Next 軸
Next 番号
For 番号 = 1 To 3 '識別番号
For 軸 = 1 To 3 '軸成分
buf01 = cls34.■除34(dt, 数字2) ' dt/2
buf01 = cls34.■加34(t, buf01) ' t+dt/2
k2vv(番号, 軸) = cls34.■乗34(dt, 速度変位成分_func(番号, 軸, buf01, r2, v2)) '速度変位_総和 m/s^2
k2rr(番号, 軸) = cls34.■乗34(dt, 位置変位成分_func(番号, 軸, buf01, r2, v2)) '位置変位成分 m/s
Next 軸
Next 番号
'k3
For 番号 = 1 To 3
For 軸 = 1 To 3
v3(番号, 軸) = cls34.■加34(v1(番号, 軸), cls34.■除34(k2vv(番号, 軸), 数字2)) 'v(識別番号,軸成分) (3,3) 前処理
r3(番号, 軸) = cls34.■加34(r1(番号, 軸), cls34.■除34(k2rr(番号, 軸), 数字2)) 'r(識別番号,軸成分) (3,3) 前処理
Next 軸
Next 番号
For 番号 = 1 To 3 '識別番号
For 軸 = 1 To 3 '軸成分
k3vv(番号, 軸) = cls34.■乗34(dt, 速度変位成分_func(番号, 軸, cls34.■加34(t, cls34.■除34(dt, 数字2)), r3, v3)) '速度変位_総和 m/s^2
k3rr(番号, 軸) = cls34.■乗34(dt, 位置変位成分_func(番号, 軸, cls34.■加34(t, cls34.■除34(dt, 数字2)), r3, v3)) '位置変位成分 m/s
Next 軸
Next 番号
'k4
For 番号 = 1 To 3
For 軸 = 1 To 3
v4(番号, 軸) = cls34.■加34(v1(番号, 軸), k3vv(番号, 軸)) 'v(識別番号,軸成分) (3,3) 前処理
r4(番号, 軸) = cls34.■加34(r1(番号, 軸), k3rr(番号, 軸)) 'r(識別番号,軸成分) (3,3) 前処理
Next 軸
Next 番号
For 番号 = 1 To 3 '識別番号
For 軸 = 1 To 3 '軸成分
k4vv(番号, 軸) = cls34.■乗34(dt, 速度変位成分_func(番号, 軸, cls34.■加34(t, dt), r4, v4)) '速度変位_総和 m/s^2
k4rr(番号, 軸) = cls34.■乗34(dt, 位置変位成分_func(番号, 軸, cls34.■加34(t, dt), r4, v4)) '位置変位成分 m/s
Next 軸
Next 番号
Dim 数字6 As String
数字6 = "6.00000000000000000000000000000000000E+0000"
Dim buf_v1 As String, buf_v1_a As String, buf_v1_b As String
Dim buf_r1 As String, buf_r1_a As String, buf_r1_b As String
For 番号 = 1 To 3
For 軸 = 1 To 3
buf_v1_a = cls34.■乗34(数字2, cls34.■加34(k3vv(番号, 軸), k2vv(番号, 軸)))
buf_v1_b = cls34.■加34(k1vv(番号, 軸), k4vv(番号, 軸))
buf_v1 = cls34.■加34(buf_v1_a, buf_v1_b)
buf_r1_a = cls34.■乗34(数字2, cls34.■加34(k2rr(番号, 軸), k3rr(番号, 軸)))
buf_r1_b = cls34.■加34(k1rr(番号, 軸), k4rr(番号, 軸))
buf_r1 = cls34.■加34(buf_r1_a, buf_r1_b)
v1(番号, 軸) = cls34.■加34(v1(番号, 軸), cls34.■除34(buf_v1, 数字6)) ' m/s
r1(番号, 軸) = cls34.■加34(r1(番号, 軸), cls34.■除34(buf_r1, 数字6)) ' m
Next 軸
Next 番号
'便宜的に戻り値を行列 ret(18)でまとめ
ret(1) = r1(1, 1)
ret(2) = r1(1, 2)
ret(3) = r1(1, 3)
ret(4) = v1(1, 1)
ret(5) = v1(1, 2)
ret(6) = v1(1, 3)
ret(7) = r1(2, 1)
ret(8) = r1(2, 2)
ret(9) = r1(2, 3)
ret(10) = v1(2, 1)
ret(11) = v1(2, 2)
ret(12) = v1(2, 3)
ret(13) = r1(3, 1)
ret(14) = r1(3, 2)
ret(15) = r1(3, 3)
ret(16) = v1(3, 1)
ret(17) = v1(3, 2)
ret(18) = v1(3, 3)
RK4 = ret
Set cls34 = Nothing
End Function
Function 速度変位成分_func(ByVal 被番号, ByVal 軸, ByVal t, ByVal R, ByVal v) As String ' （速度変位 m/s^2）を返す
Dim cls34 As clstaby_preci34
Set cls34 = New clstaby_preci34 ' インスタンスを作成
Dim sum_a As String
Dim rrr As String, rrr3 As String
Dim 差異(3) As String
Dim 基準天体 As Integer
Dim buf As String, buf2 As String
Dim 数字_1 As String
数字_1 = "-1.00000000000000000000000000000000000E+0000"
rrr = "0.00000000000000000000000000000000000E+0000" ' 初期化
sum_a = "0.00000000000000000000000000000000000E+0000" ' 初期化
' dv/dt = -Ｇ・Ｍ・(→r)／r^3 ・・・ 万有引力の法則
For 基準天体 = 1 To 3 '番号
If 被番号 <> 基準天体 Then ' 相違する場合だけ計算する
差異(1) = cls34.■減34(R(被番号, 1), R(基準天体, 1)) 'r(識別番号,軸成分) x成分
差異(2) = cls34.■減34(R(被番号, 2), R(基準天体, 2)) 'r(識別番号,軸成分) y成分
差異(3) = cls34.■減34(R(被番号, 3), R(基準天体, 3)) 'r(識別番号,軸成分) z成分
'r = √(x^2+y^2+z^2)
buf = cls34.■加34(cls34.■乗34(差異(1), 差異(1)), cls34.■乗34(差異(2), 差異(2)))
buf = cls34.■加34(buf, cls34.■乗34(差異(3), 差異(3)))
rrr = cls34.★Sqr34(buf)
' sum_a = -G * m(基準天体) * 差異(軸) / rrr / rrr / rrr + sum_a '相違し、且つ、計算する軸成分だけを積算する
buf2 = cls34.■乗34(cls34.■乗34(cls34.■乗34(数字_1, G), m(基準天体)), 差異(軸))
rrr3 = cls34.■累乗34(rrr, -3)
sum_a = cls34.■加34(cls34.■乗34(buf2, rrr3), sum_a)
End If
Next 基準天体
速度変位成分_func = sum_a '加速度 １成分のみ
Set cls34 = Nothing
End Function
Function 位置変位成分_func(ByVal 被番号, ByVal 軸, ByVal t, ByVal R, ByVal v) As String ' （位置変位 m/s）を返す
' dx/dt = v
位置変位成分_func = v(被番号, 軸) '単に速度ベクトルを返す １成分のみ
End Function