３体問題（番外編）

ピタゴラスの３体問題・・エクセルで四倍精度計算
　

-自作４倍精度クラスでの処理時間-　　　　　　　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

Math TOP