importmatplotlib.pyplotaspltimportnumpyasnpfromnumpyimportcos,sinimportmatplotlib.animationasanimationG=9.8# acceleration due to gravity, in m/s^2L1=1.0# length of pendulum 1 in mL2=1.0# length of pendulum 2 in mL=L1+L2# maximal length of the combined pendulumM1=1.0# mass of pendulum 1 in kgM2=1.0# mass of pendulum 2 in kgt_stop=2.5# how many seconds to simulatehistory_len=500# how many trajectory points to displaydefderivs(t,state):dydx=np.zeros_like(state)dydx[0]=state[1]delta=state[2]-state[0]den1=(M1+M2)*L1-M2*L1*cos(delta)*cos(delta)dydx[1]=((M2*L1*state[1]*state[1]*sin(delta)*cos(delta)+M2*G*sin(state[2])*cos(delta)+M2*L2*state[3]*state[3]*sin(delta)-(M1+M2)*G*sin(state[0]))/den1)dydx[2]=state[3]den2=(L2/L1)*den1dydx[3]=((-M2*L2*state[3]*state[3]*sin(delta)*cos(delta)+(M1+M2)*G*sin(state[0])*cos(delta)-(M1+M2)*L1*state[1]*state[1]*sin(delta)-(M1+M2)*G*sin(state[2]))/den2)returndydx# create a time array from 0..t_stop sampled at 0.02 second stepsdt=0.01t=np.arange(0,t_stop,dt)# th1 and th2 are the initial angles (degrees)# w10 and w20 are the initial angular velocities (degrees per second)th1=120.0w1=0.0th2=-10.0w2=0.0# initial statestate=np.radians([th1,w1,th2,w2])# integrate the ODE using Euler's methody=np.empty((len(t),4))y[0]=stateforiinrange(1,len(t)):y[i]=y[i-1]+derivs(t[i-1],y[i-1])*dt# A more accurate estimate could be obtained e.g. using scipy:##   y = scipy.integrate.solve_ivp(derivs, t[[0, -1]], state, t_eval=t).y.Tx1=L1*sin(y[:,0])y1=-L1*cos(y[:,0])x2=L2*sin(y[:,2])+x1y2=-L2*cos(y[:,2])+y1fig=plt.figure(figsize=(5,4))ax=fig.add_subplot(autoscale_on=False,xlim=(-L,L),ylim=(-L,1.))ax.set_aspect('equal')ax.grid()line,=ax.plot([],[],'o-',lw=2)trace,=ax.plot([],[],'.-',lw=1,ms=2)time_template='time =%.1fs'time_text=ax.text(0.05,0.9,'',transform=ax.transAxes)defanimate(i):thisx=[0,x1[i],x2[i]]thisy=[0,y1[i],y2[i]]history_x=x2[:i]history_y=y2[:i]line.set_data(thisx,thisy)trace.set_data(history_x,history_y)time_text.set_text(time_template%(i*dt))returnline,trace,time_textani=animation.FuncAnimation(fig,animate,len(y),interval=dt*1000,blit=True)plt.show()

DownloadJupyternotebook:double_pendulum.ipynb

DownloadPythonsourcecode:double_pendulum.py

Downloadzipped:double_pendulum.zip