XLPack for Python および XLPack for Matplotlib の使用例 (2) アニメーション (単振り子)

下図のように, 長さ \(l\) の糸 (重さは無視できるものとします) で吊り下げられた重さ \(m\) の重りが左右に振れているものとします. 糸が垂直な状態からの角度を \(\theta\) で表します. これを単振り子といいます.

この運動は, 時間を \(t\) として次の微分方程式で表すことができます. なお, \(g\) は重力加速度 \((9.81 m/s^2)\) です.

\(d^2\theta/dt^2 = -(mg/l)sin\theta\)

ここでは簡単のために \(m = 1\) とし, 初期条件を \(\theta = \pi/4, d\theta/dt = 0\) として \(t = 0 \sim 10\) における数値解を求め, アニメーションで表示することにします.

1. 全体の構成

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

Sub Main()
    Dim N As Long, T As Double, Y() As Double, Yp() As Double
    Dim Nt As Long, Tend As Double
    Dim X1() As Double, Y1() As Double
    Dim Info As Long
    N = 1
    ReDim Y(N - 1), Yp(N - 1)
    '-- Initial values
    T = 0
    Y(0) = Atn(1)   ' Pi() / 4
    Yp(0) = 0
    '-- Solve ODE
    Tend = 10
    Nt = 200
    ReDim X1(Nt - 1), Y1(Nt - 1)
    Call ODE_2(N, T, Y(), Yp(), Nt, Tend, X1(), Y1(), Info)
    If Info <> 0 Then Exit Sub
    '-- Display animation
    Call Animate(Nt, X1(), Y1(), Tend / Nt)
End Sub

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

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

2. 微分方程式の解

単振り子の微分方程式は 2 階常微分方程式なので, 通常は次のように 1 階の連立常微分方程式の形にして解きます.

\(d\theta_1/dt = \theta_2\)
\(d\theta_2/dt = -(mg/l)sin\theta_1\)

しかし, この場合は 1 階の項に依存しない特別な形をしているため XLPack ライブラリのサブルーチン Dopn43() (ルンゲ・クッタ・ニュストレム法) を使ってそのままの形で解くことができます.

Const L = 1, M = 1, G = 9.81
Const Mode = 3

Sub F_2(N As Long, T As Double, Y() As Double, Ypp() As Double)
    Ypp(0) = -(M * G / L) * Sin(Y(0))
End Sub

Sub ODE_2(N As Long, T As Double, Y() As Double, Yp() As Double, Nt As Long, Tend As Double, X1() As Double, Y1() As Double, Info As Long)
    Dim RTol(0) As Double, ATol(0) As Double
    Dim Dt As Double, Tout As Double, I As Long
    RTol(0) = 0.0000000001 '1e-10
    ATol(0) = RTol(0)
    Dt = Tend / Nt
    Info = 0
    For I = 0 To Nt - 1
        Tout = (I + 1) * Dt
        Call Dopn43(N, AddressOf F_2, T, Y(), Yp(), Tout, Tend, RTol(), ATol(), Mode, Info)
        If Info < 0 Or Info > 10 Then
            MsgBox ("Error in Dopn43: Info =" + Str(Info))
            Exit Sub
        End If
        X1(I) = L * Sin(Y(0))
        Y1(I) = -L * Cos(Y(0))
    Next
End Sub

得られた解は, XY 座標に変換して X1(), Y1() に返します.

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

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

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

Sub Animate(Nt As Long, X1() As Double, Y1() As Double, Dt As Double)
    Dim I As Long
    For I = 0 To Nt - 1
        Cells(6 + I, 2) = X1(I)
        Cells(6 + I, 3) = Y1(I)
    Next
    '-- Python codes
    PY "import matplotlib.pyplot as plt"
    PY "import matplotlib.animation as animation"
    PY "from XLPackPy import *"
    PY "nt =" + Str(Nt)
    PY "dt =" + Str(Dt)
    PY "x1 = xlc_get(5, 1, 4 + nt)"
    PY "y1 = xlc_get(5, 2, 4 + nt)"
    PY "fig = plt.figure()"
    PY "ax = fig.gca()"
    PY "ax.set_xlim(-1.0, 1.0)"
    PY "ax.set_ylim(-1.5, 0.5)"
    PY "ax.set_aspect('equal')"
    PY "ax.grid()"
    PY "ax.set_title('Animation (Pendulum)')"
    PY "xt = []"
    PY "yt = []"
    PY "artist = []"
    PY _
      "for i in range(nt):" + vbNewLine + _
      "  pend, = ax.plot([0, x1[i]], [0, y1[i]], 'o-', color = 'b', lw=2)" + vbNewLine + _
      "  xt.append(x1[i])" + vbNewLine + _
      "  yt.append(y1[i])" + vbNewLine + _
      "  trace, = ax.plot(xt, yt, ':c', lw=1)" + vbNewLine + _
      "  artist.append((pend, trace))"
    PY "ani = animation.ArtistAnimation(fig, artist, interval = dt * 1000)"
    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, X1() As Double, Y1() As Double, Dt As Double)
    Dim Fig As Figure, Ax As Axs, Artist() As PyObject, Anim As Animation
    Dim X(1) As Double, Y(1) As Double
    Dim I As Long
    ReDim Artist(Nt - 1, 1)
    Set Fig = Plt.Figure()
    Set Ax = Fig.Gca()
    Call Ax.Set_xlim(-1, 1)
    Call Ax.Set_ylim(-1.5, 0.5)
    Call Ax.Set_aspect(1)
    Call Ax.Grid
    Call Ax.Set_title("Animation (Pendulum)")
    X(0) = 0: Y(0) = 0
    For I = 0 To Nt - 1
        X(1) = X1(I): Y(1) = Y1(I)
        Set Artist(I, 0) = Ax.Plot(2, X(), Y(), "o-", "color=b,lw=2")
        Set Artist(I, 1) = Ax.Plot(I + 1, X1(), Y1(), ":c", "lw=1")
    Next
    Set Anim = Animation.ArtistAnimation(Fig, Nt, 2, Artist(), "interval =" + Str(Dt * 1000))
    Call Plt.Show
End Sub
公開日: 2026/04/28 (更新日: 2026/04/28)