シンプレクティック8次公式のExcel VBAのプログラムソース

一周年記念のソース公開シリーズ第二段。

まず係数は吉田春夫さんの論文　"Construction of higher order symplectic integrators", Physics Letters A,Volume 150, number 5,6,7(1990)から引っ張ってくることにした。いろいろ書かれている中で、8次の公式のSolution Cを使うことにした。論文ではwiで書いてあるのを使いやすいようにciとdiに書き直す。

i	ci	di
0	0	0
1	0.314515325	0.62903065
2	0.999190057	1.369349464
3	0.152381158	-1.064587148
4	0.299385476	1.6633581
5	-0.007805591	-1.678969283
6	-1.61921866	-1.559468038
7	-0.623838613	0.311790812
8	0.985390848	1.658990885
9	0.985390848	0.311790812
10	-0.623838613	-1.559468038
11	-1.61921866	-1.678969283
12	-0.007805591	1.6633581
13	0.299385476	-1.064587148
14	0.152381158	1.369349464
15	0.999190057	0.62903065
16	0.314515325	0

（間違っていても責任はとりません。。。）ここで、

q_i=q_i-1 + dt c_i p_i-1/m

p_i=p_i-1 - dt d_i ∂U(q_i)/∂q

プログラムはこちら。以前に書いた8の字を描く三体問題の例です。こんな感じで。

「symplectic8.txt」をダウンロード

適当に配列とかnの値を変えてもらえば。

Option Explicit

Dim qx(5) As Double, px(5) As Double
Dim qy(5) As Double, py(5) As Double
Dim n As Integer

Private Sub CommandButton1_Click()
    Dim w(7) As Double
    Dim c(16) As Double
    Dim d(16) As Double
    Dim dt As Double
    Dim i As Integer, k As Integer, j As Integer
    Dim pi As Double
    Dim ang As Double
   
    Application.ScreenUpdating = False
   
   
    w(1) = 0.311790812418427
    w(2) = -1.55946803821447
    w(3) = -1.6789692825964
    w(4) = 1.66335809963315
    w(5) = -1.06458714789183
    w(6) = 1.36934946416871
    w(7) = 0.629030650210433
    w(0) = 1# - 2# * (w(1) + w(2) + w(3) + w(4) + w(5) + w(6) + w(7))
   
    d(1) = w(7)
    d(15) = w(7)
    d(2) = w(6)
    d(14) = w(6)
    d(3) = w(5)
    d(13) = w(5)
    d(4) = w(4)
    d(12) = w(4)
    d(5) = w(3)
    d(11) = w(3)
    d(6) = w(2)
    d(10) = w(2)
    d(7) = w(1)
    d(9) = w(1)
    d(8) = w(0)
    d(16) = 0#
   
    c(1) = 0.5 * w(7)
    c(16) = 0.5 * w(7)
    c(2) = 0.5 * (w(7) + w(6))
    c(15) = 0.5 * (w(7) + w(6))
    c(3) = 0.5 * (w(6) + w(5))
    c(14) = 0.5 * (w(6) + w(5))
    c(4) = 0.5 * (w(5) + w(4))
    c(13) = 0.5 * (w(5) + w(4))
    c(5) = 0.5 * (w(4) + w(3))
    c(12) = 0.5 * (w(4) + w(3))
    c(6) = 0.5 * (w(3) + w(2))
    c(11) = 0.5 * (w(3) + w(2))
    c(7) = 0.5 * (w(2) + w(1))
    c(10) = 0.5 * (w(2) + w(1))
    c(8) = 0.5 * (w(1) + w(0))
    c(9) = 0.5 * (w(1) + w(0))
   
    n = 3
   

    qx(1) = 0.97000436
    qy(1) = -0.24308753
   
    px(1) = -0.5 * 0.93240737
    py(1) = -0.5 * 0.86473146
   
    qx(2) = -0.97000436
    qy(2) = 0.24308753
   
    px(2) = -0.5 * 0.93240737
    py(2) = -0.5 * 0.86473146
   
    qx(3) = 0#
    qy(3) = 0#
   
    px(3) = 0.93240737
    py(3) = 0.86473146
   
    Worksheets("Sheet1").Cells(1, 2) = 0#
    For j = 1 To n
        Worksheets("Sheet1").Cells(1, 1 + 2 * j) = qx(j)
        Worksheets("Sheet1").Cells(1, 2 + 2 * j) = qy(j)
    Next j

   
    dt = 2.1 / 1000
   
    For i = 1 To 1000
        Worksheets("Sheet1").Cells(i + 1, 2) = dt * i
        For k = 1 To 16
            For j = 1 To n
                qx(j) = qx(j) + dt * c(k) * px(j)
                qy(j) = qy(j) + dt * c(k) * py(j)
            Next j
            For j = 1 To n
                px(j) = px(j) + dt * d(k) * fx(j)
                py(j) = py(j) + dt * d(k) * fy(j)
            Next j
        Next k
        For j = 1 To n
            Worksheets("Sheet1").Cells(1 + i, 1 + 2 * j) = qx(j)
            Worksheets("Sheet1").Cells(1 + i, 2 + 2 * j) = qy(j)
        Next j

    Next i

    Application.ScreenUpdating = True
   
End Sub

Function fx(k As Integer) As Double
    Dim i As Integer
    Dim f As Double, R As Double
    Dim G As Double
    G = 1
   
    f = 0
    For i = 1 To n
        If i <> k Then
            R = Sqr((qx(i) - qx(k)) ^ 2 + (qy(i) - qy(k)) ^ 2)
            f = f + G * (qx(i) - qx(k)) / (R ^ 3)
        End If
    Next i
    fx = f

End Function

Function fy(k As Integer) As Double
    Dim i As Integer
    Dim f As Double, R As Double
    Dim G As Double
    G = 1
      
    f = 0
    For i = 1 To n
        If i <> k Then
            R = Sqr((qx(i) - qx(k)) ^ 2 + (qy(i) - qy(k)) ^ 2)
            f = f + G * (qy(i) - qy(k)) / (R ^ 3)
        End If
    Next i
    fy = f

End Function