« 売布神社へ行ってきた。 | トップページ | ソフトバンクのCMでさだまさしが歌う曲の歌詞を全部読むと実は結構いい。 »

2009年4月18日 (土)

シンプレクティック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

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

qi=qi-1 + dt ci pi-1/m

pi=pi-1 - dt di ∂U(qi)/∂q

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

8noji3body

「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

« 売布神社へ行ってきた。 | トップページ | ソフトバンクのCMでさだまさしが歌う曲の歌詞を全部読むと実は結構いい。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/44650071

この記事へのトラックバック一覧です: シンプレクティック8次公式のExcel VBAのプログラムソース:

« 売布神社へ行ってきた。 | トップページ | ソフトバンクのCMでさだまさしが歌う曲の歌詞を全部読むと実は結構いい。 »

最近のコメント

2018年11月
        1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30  
フォト
無料ブログはココログ