这里是how cart pendulum looks like
假设你有一个4个微分方程，它代表了一个动力系统的运动（推车上的钟摆），你用 scipy.integrate.odeint 解了这些方程10秒，间隔为0.01秒。
最后，你得到了大小为（1000，4）的解矩阵。对于每个差分方程，您可以获得1000个数据点。到目前为止一切都很好。例如，如果我绘制的议案之一，我可以得到美丽的图形.（* 下面的图像显示了摆杆的运动（振荡）*）
这里是Graph of theta angle
但是，而不是无聊的图形，我想做一个动画，显示了运动的车史蒂夫布伦顿没有它作为下面的链接与使用Matlab。这里是link of the cart-pend video！

为了使这些图形动画化，我实际上试图用Python做Steve Brunton在Matlab中做的事情。但结果只是一个冻结的数字，而不是移动一个。实际上，如果我从Spyder IDE运行这个脚本，我会在IPython控制台中得到1000个数字。（* 每个数字代表系统瞬时运动的快照，这很好。但我只想要一个有1000个连续帧的数字。*）
这是snap of frozen cart-pend
我写了两个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)

型