XLPack for Python および XLPack for Matplotlib の使用例 (3) アニメーション (アレンストーフ軌道)

質量が \(\mu\) と \(\mu’ (= 1 – \mu)\) の2つの天体が平面内回転運動をしており, その2つに比べて質量が無視できる第3の物体が同一平面内を運動するものとします. その方程式は次のようになります (参考文献: E. Hairer, 他「常微分方程式の数値解法 I」(2007) シュプリンガー・ジャパン).
\[
\begin{align}
& y_1” = y_1 + 2y_2′ – \mu'(y_1 + \mu)/D_1 – \mu(y_1 – \mu’)/D_2 \\
& y_2” = y_2 – 2y_1′ – \mu’y_2/D_1 – \mu y_2/D_2 \\
& ただし, \\
& D_1 = ((y_1 + \mu)^2 + y_2^2)^{3/2} \\
& D_2 = ((y_1 – \mu’)^2 + y_2^2)^{3/2} \\
& \mu = 0.012277471, \mu’ = 1 – \mu \\
\end{align}
\] 次の初期値を用いると, \(t = 17.0652165601579625588917206249\) で周期解となることが知られています.
\[
\begin{align}
& y_1(0) = 0.994, y_1′(0) = 0 \\
& y_2(0) = 0, y_2′(0) = -2.00158510637908252240537862224 \\
\end{align}
\] これはアレンストーフ軌道とよばれています.

ここでは, 1周期分の数値解を求めアニメーションで表示することにします.

1. 全体の構成

メインプログラムは以下のとおりです.

Const Tend = 17.065216560158

Sub Main()
    Dim Xt() As Double, Yt() As Double
    Dim Dt As Double, Nt As Long, Info As Long
    Dt = 0.1
    Nt = Int(Tend / Dt) + 1
    ReDim Xt(Nt), Yt(Nt)
    Call ODE(Tend, Dt, Nt, Xt(), Yt(), Info)
    Call Animate(Nt, Xt(), Yt(), Dt)
End Sub

サブルーチン ODE() は微分方程式の解を求めます. 得られた解はサブルーチン Animate() により表示されます. Animate() では XLPack for Python または XLPack for Matplotlib を使用してアニメーションで表示を行います.

実行結果は次のようになります.

2. 微分方程式の解

この微分方程式は 2 階常微分方程式なので, 1 階の連立常微分方程式の形に変数変換して 4 元の方程式として解きます.

Const Xinit = 0.994

Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
    Const Mu = 0.012277471, Mup = 1 - Mu
    Dim D1 As Double, D2 As Double
    Yp(0) = Y(2)
    Yp(1) = Y(3)
    D1 = ((Y(0) + Mu) ^ 2 + Y(1) ^ 2) ^ (3# / 2#)
    D2 = ((Y(0) - Mup) ^ 2 + Y(1) ^ 2) ^ (3# / 2#)
    Yp(2) = Y(0) + 2 * Y(3) - Mup * (Y(0) + Mu) / D1 - Mu * (Y(0) - Mup) / D2
    Yp(3) = Y(1) - 2 * Y(2) - Mup * Y(1) / D1 - Mu * Y(1) / D2
End Sub

Sub ODE(Tend As Double, Dt As Double, Nt As Long, Xt() As Double, Yt() As Double, Info As Long)
    Const N = 4, Mode = 3, Tol = 0.0000000001 '1.0e-10
    Dim Y(N - 1) As Double, T As Double, Tout As Double
    Dim RTol(0) As Double, ATol(0) As Double, I As Integer
    RTol(0) = Tol
    ATol(0) = Tol
    T = 0
    Y(0) = Xinit: Y(1) = 0: Y(2) = 0: Y(3) = -2.00158510637908
    Info = 0
    Xt(0) = Y(0)
    Yt(0) = Y(1)
    For I = 0 To Nt - 1
        Tout = (I + 1) * Dt
        If Tout > Tend Then Tout = Tend
        Call Derkfa(N, AddressOf F, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info)
        If Info < 0 Or Info > 10 Then Exit Sub
        Xt(I + 1) = Y(0)
        Yt(I + 1) = Y(1)
    Next
End Sub

得られた解は, Xt(), Yt() に保存します.

3. 解のアニメーション表示

得られた解を 表示する Python のコードを XLPack for Python を使って直接記述することができます.

配列データはそのまま渡すことができないので, いったんワークシートに書き出してから Python 側で xlc_get() を使って読み込んでいます.

Sub Animate(Nt As Long, Xt() As Double, Yt() As Double, Dt As Double)
    Const Col = 2, Row = 6
    Dim I As Long
    For I = 0 To Nt
        Cells(Row + I, Col) = Xt(I)
        Cells(Row + I, Col + 1) = Yt(I)
    Next
    '-- Python codes
    PY "import matplotlib.pyplot as plt"
    PY "import matplotlib.animation as animation"
    PY "from XLPackPy import *"
    PY "fig = plt.figure()"
    PY "ax = fig.gca()"
    PY "ax.grid()"
    PY "ax.plot(0, 0, '-o', color = 'orange')"
    PY "ax.plot(1, 0, '-o', color = 'yellow')"
    PY "ax.set_xlim(-1.5, 1.5)"
    PY "ax.set_ylim(-1.5, 1.5)"
    PY "xt = xlc_get(5, 1," + Str(5 + Nt) + ")"
    PY "yt = xlc_get(5, 2," + Str(5 + Nt) + ")"
    PY "xtrace = []; ytrace = []; artist = []"
    PY "for i in range(len(xt)):" + vbNewLine + _
        "  a_x = [xt[i]]; a_y = [yt[i]]" + vbNewLine + _
        "  a, = ax.plot(a_x, a_y, 'o-', color = 'b')" + vbNewLine + _
        "  xtrace.append(xt[i]); ytrace.append(yt[i])" + vbNewLine + _
        "  trace, = ax.plot(xtrace, ytrace, ':c', lw=1)" + vbNewLine + _
        "  artist.append((a, trace))"
    PY "ani = animation.ArtistAnimation(fig, artist, interval = 100, blit = True, repeat = False)"
    PY "plt.show()"
End Sub

XLPack for Matplotlib を使うと同じ処理を全て VBA で記述することができます. 内部的には Python を呼び出して処理します.

Dim Plt As New PyPlot
Dim Animation As New Animation

Sub Animate(Nt As Long, Xt() As Double, Yt() As Double, Dt As Double)
    Dim Ex(0) As Double, Ey(0) As Double, Mx(0) As Double, My(0) As Double
    Dim X(0) As Double, Y(0) As Double, I As Long
    Dim Fig As Figure, Ax As Axs, Artist() As PyObject
    Dim Emarker As PyObject, Mmarker As PyObject
    Dim Anim As ArtistAnimation
    Ex(0) = 0: Ey(0) = 0
    Mx(0) = Xinit: My(0) = 0
    ReDim Artist(Nt, 1)
    Set Fig = Plt.Figure()
    Set Ax = Fig.Gca()
    Call Ax.Grid
    Call Ax.Set_xlim(-1.5, 1.5)
    Call Ax.Set_ylim(-1.5, 1.5)
    Set Emarker = Ax.Plot(1, Ex(), Ey(), "", "color=orange,marker=o")
    Set Mmarker = Ax.Plot(1, Mx(), My(), "", "color=yellow,marker=o")
    For I = 0 To Nt
        X(0) = Xt(I): Y(0) = Yt(I)
        Set Artist(I, 0) = Ax.Plot(1, X(), Y(), "", "color=b,marker=o")
        Set Artist(I, 1) = Ax.Plot(I + 1, Xt(), Yt(), ":c", "lw=1")
    Next
    Dim Ani As PyObject
    Set Anim = Animation.ArtistAnimation(Fig, Nt + 1, 2, Artist(), "interval =" + Str(Dt * 1000) + ", repeat = False")
    Call Plt.Show
End Sub
公開日: 2026/04/28 (更新日: 2026/04/28)