How to animate motion of a dynamical system in Python?

Question

Here's the how cart pendulum looks like

Imagine you have a 4 differantial equations which represents the motion of a dynamic system (pendulum on a cart) and you solved these equations using scipy.integrate.odeint for 10 seconds with interval of 0.01 seconds.

Finally you get the solution matrix with size (1000,4). For each diff eqs you get 1000 data points. Everything is ok so far. For example, If I plot one of the motion I can get beautiful graphics.(Below image shows the motion of the pendulum rod(oscillating))

Here's Graph of theta angle

But, instead of boring graphics and I want to make an animation that shows the motion of the cart as Steve Brunton did it as below link with using Matlab. Here's link of the cart-pend video!

====================================================================

To animate the figures I actually tried to do what Steve Brunton did in Matlab, with Python. But the result is just a frozen figure instead of moving one. Actually If I run this script from Spyder IDE, I get 1000 figures in the IPython console.(Each figure represents a snapshot of the system's instantaneous motion which is good. But I want just one figure with 1000 sequantial frames on it.)

Here's the snap of frozen cart-pend

I've written two python scripts. One for only plotting the other is for solving the diff eqs and feed the results to the other one.

~~~~~~~~~~~~~~~~~~~~~~~~~

This code is for plotting the animated figures.

from math import sqrt, sin, cos
import matplotlib.pyplot as plt
from matplotlib import animation

def draw_cart(states, m, M, L):

    x = states[0]       # Position of the center of the cart
    theta = states[3]   # Angle of the pendulum rod

    #Dimensions
    W = 1*sqrt(M/5)     # Cart width
    H = .5*sqrt(M/5)    # Cart Height
    wr = .2             # Wheel radius
    mr = .3*sqrt(m)     # Mass Radius

    #Positions

    y = wr/2+ H/2       # Cart Vertical Position

    w1x = x-.9*W/2      # Left Wheel x coordinate
    w1y = 0             # Left wheel y coordinate
    w2x = x+(.9*W/2)    # Right Wheel x coordinate
    w2y = 0             # Right Wheel y coordinate

    # Pendulum Mass x-y coordinates
    px = x+(L*sin(theta))
    py = y-(L*cos(theta))



    #Identfying Figure
    plt.figure()
    plt.axes(xlim=(-5, 5), ylim=(-2, 2.5))

    # Plotting the base line
    line = plt.Line2D((-10, 10), (0, 0), color='k', linewidth=2)
    plt.gca().add_line(line)
    plt.hold(True)


    # Shapes
    rectangle1 = plt.Rectangle((x-(W/2), (y-H/2)), W, H, fill=True, color='b') # Cart

    rectangle2= plt.Rectangle((px-(mr/2), py-(mr/2)), mr, mr, fill=True, color='r') # Pendulum mass

    circle2 = plt.Circle((w1x, w1y), wr/2, fill=True, color='g') #Left whell

    circle3 = plt.Circle((w2x, w2y), wr/2, fill=True, color='g')  #Right whell

    plt.plot((x, px), (y, py), 'k', lw=2) #Pendulum rod

    #Adding shapes to the figure
    plt.gca().add_patch(rectangle1)
    plt.gca().add_patch(rectangle2)
    plt.gca().add_patch(circle2) 
    plt.gca().add_patch(circle3) 

    # Showing the figure
    plt.show()

    plt.hold(False)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This is the other code for solving the diff eqs and feeding the solution to the above code.

from math import pi, sin, cos
import numpy as np
from scipy.integrate import odeint
import draw_cart_pend_rt
import matplotlib.pyplot as plt

# System Parameters
m = 1
M = 5
L = 2
g = -10
d = 1
u = 0


def cart_pend_dynamics(states, tspan):

    Sy = sin(states[2])
    Cy = cos(states[2])
    D = m*L*L*(M+(m*(1-(Cy**2))))

    state_derivatives = np.zeros_like(states)

    state_derivatives[0] = states[1]
    state_derivatives[1] = ((1/D)*(((-m**2)*(L**2)*g*Cy*Sy)+(m*(L**2)*(m*L*(states[3]**2)*Sy-d*(states[1])))))+(m*L*L*(1/D)*u)
    state_derivatives[2] = states[3]
    state_derivatives[3] = ((1/D)*((m+M)*m*g*L*Sy-m*L*Cy*(m*L*(states[3])**2*Sy-d*states[1])))-(m*L*Cy*(1/D)*u)+(0.01*1)

    return state_derivatives



def solution_of_cartpend(dt):

    # Initial conditions to solve diff eqs
    states = np.array([0.0, 0.0, pi, 0.5])  # Left to right, cart; position-velocity, pend mass; angle-angular velocity

    tspan = np.arange(0, 10, dt)

    state_sol = odeint(cart_pend_dynamics, states, tspan)

    return state_sol


# Time Interval
dt = 0.01

solution = solution_of_cartpend(dt)

x_den, y_den = solution.shape



# Validating the solution
plt.axes(xlim=(0,10), ylim=(-10,10))
t = np.arange(0, 10, dt)
plt.gca().plot(t, (solution[:, 2]), 'b', label='theta1')

# Animating the figures
for i in range(x_den):

    draw_cart_pend_rt.draw_cart(solution[i,:], m, M, L)

Answer 1

Please try the following code.

import numpy
import scipy
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
#from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation


def cart_pend(t, x, pend_mass, cart_Mass, L, g, d, u):
    """

    This function is to idenstify the our equations that link our states to the inputs

    """

    Sx = numpy.sin(x[2])
    Cx = numpy.cos(x[2])
    D = pend_mass * L * L * (cart_Mass + pend_mass * (1 - Cx**2))

    dx1_dt = x[1]
    dx2_dt = (1 / D) * (
        -(pend_mass**2) * L**2 * g * Cx * Sx
        + pend_mass * L**2 * (pend_mass * L * x[3] ** 2 * Sx - d * x[1])
    ) + pend_mass * L * L * (1 / D) * u
    dx3_dt = x[3]
    dx4_dt = (1 / D) * (
        (pend_mass + cart_Mass) * pend_mass * g * L * Sx
        - pend_mass * L * Cx * (pend_mass * L * x[3] ** 2 * Sx - d * x[1])
    ) - pend_mass * L * Cx * (1 / D) * u

    return [dx1_dt, dx2_dt, dx3_dt, dx4_dt]


m = 1
M = 5
L = 2
g = -10
d = 1
tspan = (1.0, 10.0)
x0 = [0, 0, numpy.pi, 0.5]
print("the shape of initial state vector is", numpy.shape(x0))
sol = solve_ivp(lambda t, x: cart_pend(t, x, m, M, L, g, d, -1), tspan, x0)
t = sol.t
y1, y2, y3, y4 = sol.y

print("Time points:", t)
print("Solution for y1:", y1)
print("Solution for y2:", y2)

plt.plot(t, y1, label="y1")
plt.plot(t, y2, label="y2")
plt.xlabel("Time")
plt.ylabel("State Variables")
plt.legend()
plt.grid(True)
plt.show()

#fig, ax = plt.subplots()
#(line,) = ax.plot([], [], "b")
#
#
#def init():
#    ax.set_xlim(-2, 2)
#    ax.set_ylim(-2, 2)
#    return (line,)
#
#
#def update(frame):
#    line.set_data(y1[:frame], y2[:frame])
#    return (line,)
#
#
#ani = FuncAnimation(fig, update, frames=len(t), init_func=init, blit=True)
#ani.save("simulation.gif", writer="Pillow")
#
#plt.show()


# Initialize figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)

# Initialize shapes
shape1 = plt.Circle((0, 0), 0.1, color='red')
shape2 = plt.Rectangle((0, 0), 0.2, 0.2, color='blue')

# Update function to compute new positions
def update(frame):
    theta = sol.y[0, frame]
    x = sol.y[2, frame]
    shape1.center = (numpy.sin(theta), numpy.cos(theta))
    shape2.set_xy([x - 0.1, -0.1])
    return shape1, shape2

# Animation function
def animate(frame):
    ax.clear()
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.add_patch(shape1)
    ax.add_patch(shape2)
    return update(frame)

# Create animation
anim = animation.FuncAnimation(fig, animate, frames=sol.y.shape[1], interval=100)

# Display the animation
plt.show()

# Convert animation to HTML format
html_anim = anim.to_jshtml()

# Save animation to an HTML file
with open("animation.html", "w") as f:
    f.write(html_anim)

# Display the animation
plt.close()