Example program of XLPack for Python and XLPack for Matplotlib: (2) Animation (Simple pendulum)

As shown in the figure below, consider a mass of weight \(m\) suspended by a string of length \(l\) (whose weight is assumed to be negligible) oscillating from side to side. Let \(\theta\) denote the angle measured from the vertical position of the string. This system is called a simple pendulum.

This motion can be described by the following differential equation, where \(t\) denotes time. Here, \(g\) is the acceleration due to gravity \((9.81 m/s^2)\).

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

For simplicity, we assume \(m = 1\), and compute the numerical solution over the interval \(t = 0 \sim 10\) with the initial conditions \(\theta = \pi/4\) and \(d\theta/dt = 0\). The result will then be displayed as an animation.

1. Overall structure

The main program is as follows.

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

The subroutine ODE_2() solves the differential equation. The obtained solution is displayed by the subroutine Animate(). In Animate(), animation is performed using XLPack for Python or XLPack for Matplotlib.

The execution result will be as follows.

2. Solutions to differential equations

Since the differential equation of a simple pendulum is a second-order ordinary differential equation, it is usually rewritten as the following system of first-order ordinary differential equations before being numerically solved.

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

However, in this case, the equation has a special form that does not depend on first-order terms, so it can be solved directly in its original form using the subroutine Dopn43() (Runge-Kutta-Nystrom method) in the XLPack library.

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

The obtained solution is converted into XY coordinates and returned in X1() and Y1().

3. Animation of the Solution

The obtained solution can be directly visualized using Python code written with using XLPack for Python.

Since array data cannot be passed directly, the data is first written to a worksheet and then read on the Python side using 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 "anim = animation.ArtistAnimation(fig, artist, interval = dt * 1000)"
    PY "plt.show()"
End Sub

Using XLPack for Matplotlib, the same processing can be written entirely in VBA. Internally, the processing is carried out by calling 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