harmonic motion animation with python

Question

I am trying to plot an animation of a simple pendulum using the model https://matplotlib.org/gallery/animation/double_pendulum_sgskip.html.

My code is as follows :

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()

EDIT : I am getting the error message :

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

  <Figure size 432x288 with 1 Axes>

I have checked what my function step_solver is printing with

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

and after recommendation the other y output for odeint which gave only 2 values for my ysol and ysolv variable for 2200 time values.

I was expecting to get a range of y values for each time step. I am not sure how to sort this. Is my function step_solver wrongly coded? Why am I only getting 2 values, how can I animate the solution in the same way it was done for the double pendulum?

Any suggestions on where the problem may lie? Many thanks in advance.

Maybe you are not using an interactive backend? How do you run the code? — ImportanceOfBeingErnest, Feb 21 '19 at 00:17
@ImportanceOfBeingErnest thank you for pointing that out.Indeed I was missing it. I am still having issues though I will edit the code accordingly now — kevint, Feb 21 '19 at 01:12
Yes, `integrate.odeint` returns two objects. You are only interested in the first one of those 2 I suppose. — ImportanceOfBeingErnest, Feb 21 '19 at 01:26
@ImportanceOfBeingErnest I am not understanding how can I then populate my animate function such that I am plotting the solutions to these ODE against time elapsed? Apologies for the confusion and thank you for your effeorts — kevint, Feb 21 '19 at 02:07
Sorry, what I mean is that `odeint` returns two objects, one of which you want to ignore. That is, you want to unpack it like `desiredstate, garbage = integrate.odeint(...)` — ImportanceOfBeingErnest, Feb 21 '19 at 02:14
@ImportanceOfBeingErnest I am still struggling I have edited the code accordingly and checked the second output which has also a size of 2. I am not sure how can I plot this in the same way it was done in the link above. How can I unpack my range of y ode solutions ? — kevint, Feb 21 '19 at 03:22
You did not follow the code in the link. Can you tell what motivated the changes in the design? And what exactly do you expect to see in the plot at any given time? — Lutz Lehmann, Feb 21 '19 at 06:48
@ImportanceOfBeingErnest : `odeint` returns a list of the sample values according to the time array that was passed. If as here the time is `[t0,tf]`, then the list contains `[y0,yf]` where `yf=y(tf)` and each of these samples is a 2 vector. The computed, and thus most likely more desired, value is the second one. — Lutz Lehmann, Feb 21 '19 at 07:37
@LutzL I wanted to animate a simple pendulum as given by the above equation hence the changes above. I expect to see a pendulum oscillating back and forth — kevint, Feb 21 '19 at 07:57
This means you want to draw a line with handles like `plot([0,l*sin(phi)],[0,-l*cos(phi)], '-o', lw=2, ms=8)` where `phi=state[i,0]`? — Lutz Lehmann, Feb 21 '19 at 08:29

Lutz Lehmann · Accepted Answer · 2019-02-21T08:36:17.363

Let's follow the code from the example more closely. Which means to use just one integration call and then using slices of the result of that call in the animation.

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()

For a variant using a step generator using yield see my answer in Error in RK4 algorithm in Python. This allows to encapsulate the data for the integration loop without defining a class for it. However it is not clear how a function graph with two samples would be helpful, even if it is animated.

To animate the pendulum itself, use

##===================================
##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()

harmonic motion animation with python

1 Answers1