python的谐波运动动画

Question

我正在尝试绘制一个简单钟摆的动画 使用模型https://matplotlib.org/gallery/animation/double_pendulum_sgskip.html。

我的代码如下：

from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation

#some constants 
g = 9.81 
l = 0.1 
m = 0.01 

def sh(r,t):
    theta = r[0]
    omega = r[1]
    sh_theta = omega
    sh_omega = -g/l*sin(theta) 
    return np.array([sh_theta,sh_omega],float)

init_state = [np.radians(89.0),0]
time = np.arange(0,50.0,0.025)
time_elapsed = time

def step_solver(eq, ist, dt):

    """
    Execute one time step of length dt and update status

    """
    global time_elapsed 

    state1,state2 = integrate.odeint(eq,ist,[0,dt])
    time_elapsed += dt

    return state1, state2,time_elapsed
dt = 1/30
ysol,ysolv,timex = step_solver(sh,init_state,dt)
print("This is the y0 values: ", ysol,"y values",ysolv,"This is the time ", timex) 

##===================================
##Setting up figure and animation 
#======================================
fig = plt.figure()
ax = plt.axes(xlim = (0,2), ylim = (-2,2))
line, = ax.plot([],[],lw=2)
#time_text = ax.text(0.02,0.95,'',transform = ax.transAxes) 
#==========================================
##initialisation function: plot the background of each frame 
def init():
    line.set_data([],[])
    #time_text.set_text('')
    return line, 

def animate(ysol,timex):

    x = timex 
    y = ysol

    line.set_data(x,y)
    #time_text.set_text(str(i))
    return line,
#================================================
#Call the animator    
#================================================
anim = animation.FuncAnimation(fig,animate(ysolv,timex), init_func = init, frames = 200, interval =20, blit = True)
anim.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()

编辑： 我收到错误消息：

  ValueError: shape mismatch: objects cannot be broadcast to a single shape

  <Figure size 432x288 with 1 Axes>

我检查了我的函数 step_solver 正在打印什么

 print(len(ysol),len(ysolv),len(timex))

在推荐 odeint 的另一个 y 输出之后，它只为我的 ysol 和 ysolv 变量提供了 2200 个时间值的 2 个值。

我希望每个时间步都能得到一个 y 值范围。 我不知道如何排序。我的函数 step_solver 是否编码错误？ 为什么我只得到 2 个值，如何以与双摆相同的方式为解决方案设置动画？

关于问题可能出在哪里的任何建议？ 提前谢谢了。

Answer 1

让我们更密切地关注示例中的代码。这意味着只使用一个集成调用，然后在动画中使用该调用结果的切片。

from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation

#some constants 
g = 9.81 
l = 0.1 
m = 0.01 

def sh(r,t):
    theta, omega = r
    sh_theta = omega
    sh_omega = -g/l*sin(theta) 
    return np.array([sh_theta,sh_omega],float)

init_state = np.radians([89.0,0])
dt = 1.0/30
time = np.arange(0,50.0,dt)

state = integrate.odeint(sh,init_state,time)

##===================================
##Setting up figure and animation 
#======================================
fig = plt.figure()
ax = plt.axes(xlim = (0,10), ylim = (-2,2))
line, = ax.plot([],[],lw=2)
#==========================================
##initialisation function: plot the background of each frame 
def init():
    return line, 

def animate(i):

    x = time[i:i+30] 
    y = state[i:i+30,0]

    line.set_data(x,y)
    return line,
#================================================
#Call the animator    
#================================================
anim = animation.FuncAnimation(fig,animate, init_func = init, frames = 200, interval =20, blit = True)
plt.show()

对于使用yield 的步骤生成器的变体，请参阅我在Error in RK4 algorithm in Python 中的回答。这允许封装集成循环的数据，而无需为其定义类。然而，即使是动画，包含两个样本的函数图是否会有帮助仍不清楚。

要为钟摆本身设置动画，请使用

##===================================
##Setting up figure and animation 
#======================================
fig = plt.figure(figsize=(8,6))
ax = plt.axes(xlim = (-2*l,2*l), ylim = (-2*l,l))
line, = ax.plot([],[],'-o',lw=2,ms=8)
#==========================================
##initialisation function: plot the background of each frame 
def init():
    return line, 

def animate(i):

    phi = state[i,0]

    line.set_data([0,l*sin(phi)],[0,-l*cos(phi)])
    return line,
#================================================
#Call the animator    
#================================================
anim = animation.FuncAnimation(fig,animate, init_func = init, frames = len(time), interval =100, blit = True)
plt.show()