ピタゴラスの3体問題・・エクセルで四倍精度計算
-自作4倍精度クラスでの処理時間- 34桁関数モジュール(clstaby_preci34)
ほぼ同一のプログラムで d=0.0001 t=1 の条件で処理時間を比較すると
約2125倍遅くなる。
- 結果 -
倍精度時
1分06秒
自作4倍精度クラス使用 35時間57分54秒
(ムーアの法則だと10年で1000倍なので11年分処理が遅くなる。 注意:改善余地あり。)
◆◆◆◆ ピタゴラス3体問題 ◆◆◆◆ ※別途クラスモジュール: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"
' 質量1
m(2) =
"4.00000000000000000000000000000000000E+0000"
' 質量2
m(3) =
"5.00000000000000000000000000000000000E+0000"
' 質量3
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))
'M1 x
.Cells(cnt1 + 10, 6).Value = cls34.■34↓15(rr(1,
2))
'M1 y
.Cells(cnt1 + 10, 7).Value = cls34.■34↓15(rr(2,
1))
'M2 x
.Cells(cnt1 + 10, 8).Value = cls34.■34↓15(rr(2,
2))
'M2 y
.Cells(cnt1 + 10, 9).Value = cls34.■34↓15(rr(3,
1))
'M3 x
.Cells(cnt1 + 10, 10).Value = cls34.■34↓15(rr(3,
2))
'M3 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 = -G・M・(→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
'加速度 1成分のみ
Set cls34 = Nothing
End Function
Function 位置変位成分_func(ByVal 被番号, ByVal 軸, ByVal t, ByVal R, ByVal
v) As String
' (位置変位 m/s)を返す
' dx/dt =
v
位置変位成分_func = v(被番号, 軸)
'単に速度ベクトルを返す 1成分のみ
End Function