Как оживить движение динамической системы в Python?

Question

Вот как выглядит маятник тележки

Представьте, что у вас есть 4 дифференциальных уравнения, которые представляют движение динамической системы ( маятник на тележке ) иВы решили эти уравнения, используя scipy.integrate.odeint в течение 10 секунд с интервалом 0,01 секунды.

Наконец, вы получаете матрицу решения с размером (1000,4).За каждое различие вы получаете 1000 точек данных.Пока все в порядке.Например, если я нарисую одно из движений, я получу красивую графику. ( На рисунке ниже показано движение маятникового стержня (колебание) )

Вот График тета-угла

Но вместо скучной графики я хочу сделать анимацию, которая показывает движение тележки, как это сделал Стив Брантон, как показано ниже, по ссылке с использованием Matlab.Вот ссылка на видео в ожидании корзины !

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

Для анимации фигур яна самом деле пытался сделать то, что сделал Стив Брантон в Matlab с Python.Но результат - просто замороженная фигура, а не движущаяся.На самом деле, если я запускаю этот сценарий из Spyder IDE, я получаю 1000 фигур в консоли IPython. ( Каждая цифра представляет собой снимок мгновенного движения системы, что хорошо. Но мне нужна только одна фигура с 1000 последовательных кадров на ней.)

Вот снимок замороженной корзины-ожидания

Я написал два скрипта Python.Один только для построения графика, другой - для решения дифференциальных уравнений и передачи результатов другому.

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

Этот код предназначен для построения анимированных фигур.

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)

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

Это другой код для решения дифференциальных уравнений и подачи решения в приведенный выше код.

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)