当我尝试在python中模拟双摆时出现问题运行代码给我错误信息

问题描述

我想在python中模拟一个双摆，我在下面编写了代码，但是现在我完全被卡住了。每次运行代码时都会出错，并且已经看了几个小时了，但是我不明白为什么会出错并且模拟不起作用

from matplotlib import animation
from pylab import *
from scipy import optimize

G = 9.8  # gravitational acceleration


# Kinetic energy
def Ekin(osc):
    return 1 / (2.0 * osc.m * osc.L * osc.L) * (
        osc.p1 * osc.p1 + 2.0 * osc.p2 * osc.p2 - 2.0 * osc.p1 * osc.p2 * cos(osc.q1 - osc.q2)) / (
               1 + (sin(osc.q1 - osc.q2)) ** 2)


# Potential energy
def Epot(osc):
    return osc.m * G * osc.L * (3 - 2 * math.cos(osc.q1) - math.cos(osc.q2))


# Class that holds the parameter and state of a double pendulum
class Oscillator:
    m = 1  # mass of the pendulum bob
    L = 1  # arm length
    t = 0  # the time
    E = 15
    q1 = 0
    q2 = np.pi/2
    p1=0
    p2 = math.sqrt((E-G*(3-2*cos(q1)-cos(q2)))*(1+(sin(q1-q2))**2)*m*L**2)


class Observables:
    def __init__(self):
        self.time = []  # list to store time
        self.q1list = []  # list to store q1
        self.q2list = []  # list to store q2
        self.p1list = []  # list to store p1
        self.p2list = []  # list to store p2
        self.epot = []  # list to store potential energy
        self.ekin = []  # list to store kinetic energy
        self.etot = []  # list to store total energy
        self.poincare_q1 = []  # list to store q1 for Poincare plot
        self.poincare_p1 = []  # list to store p1 for Poincare plot


# Derivate of H with respect to p1
def dHdp1(q1,q2,p1,p2,m,L):
    return (p1-p2)*cos(q1-q2)/(m*L**2*(1+(sin(q1-q2))**2))

# Derivate of H with respect to p2
def dHdp2(q1,L):
    return (2*p2-p1*cos(q1-q2))/(m*L**2*(1+(sin(q1-q2))**2))

# Derivate of H with respect to q1
def dHdq1(q1,L):
    return 1 / (2.0 * m * L * L) * (
        -2 * (p1 * p1 + 2 * p2 * p2) * cos(q1 - q2) + p1 * p2 * (4 + 2 * (cos(q1 - q2)) ** 2)) * sin(
    q1 - q2) / (1 + (sin(q1 - q2)) ** 2) ** 2 + m * G * L * 2.0 * sin(q1)


# Derivate of H with respect to q2
def dHdq2(q1,L):
    return 1 / (2.0 * m * L * L) * (
        2 * (p1 * p1 + 2 * p2 * p2) * cos(q1 - q2) - p1 * p2 * (4 + 2 * (cos(q1 - q2)) ** 2)) * sin(q1 - q2) / (
               1 + (sin(q1 - q2)) ** 2) ** 2 + m * G * L * sin(q2)


class BaseIntegrator:
    dt = 0.01  # time step

    def integrate(self,osc,obs,):
    """ Perform a single integration step """
        self.timestep(osc,obs)

    """ Append observables to their lists """
        obs.time.append(osc.t)
        obs.q1list.append(osc.q1)
        obs.q2list.append(osc.q2)
        obs.p1list.append(osc.p1)
        obs.p2list.append(osc.p2)
        obs.epot.append(Epot(osc))
        obs.ekin.append(Ekin(osc))
        obs.etot.append(Epot(osc) + Ekin(osc))

    def timestep(self,obs):
    """ Implemented by the child classes """
        pass




# Runge-Kutta 4 integrator
class RK4Integrator(BaseIntegrator):
    def timestep(self,obs):
        dt = self.dt
        dp1dt = -dHdq1(osc.q1,osc.q2,osc.p1,osc.p2,osc.m,osc.L); dp2dt = -dHdq2(osc.q1,osc.L); dq1dt = dHdp1(osc.q1,osc.L); dq2dt = dHdp2(osc.q1,osc.L)
        k1p1 = dt*dp1dt; k1p2 = dt*dp2dt; k1q1 = dq1dt * dt; k1q2 = dq2dt * dt
        k2p1 = dt * (-dHdq1(osc.q1+k1q1/2,osc.q2+k1q2/2,osc.p1+k1p1/2,osc.p2+k1p2/2,osc.L))
        k2p2 = dt * (-dHdq2(osc.q1+k1q1/2,osc.L))
        k2q1 = dt * dHdp1(osc.q1+k1q1/2,osc.L)
        k2q2 = dt * dHdp2(osc.q1+k1q1/2,osc.L)
        k3p1 = dt * (-dHdq1(osc.q1+k2q1/2,osc.q2+k2q2/2,osc.p1+k2p1/2,osc.p2+k2p2/2,osc.L))
        k3p2 = dt * (-dHdq2(osc.q1+k2q1/2,osc.L))
        k3q1 = dt * dHdp1(osc.q1+k2q1/2,osc.L)
        k3q2 = dt * dHdp2(osc.q1+k2q1/2,osc.L)
        k4p1 = dt * (-dHdq1(osc.q1+k3q1,osc.q2+k3q2,osc.p1+k3p1,osc.p2+k3p2,osc.L))
        k4p2 = dt * (-dHdq2(osc.q1+k3q1,osc.L))
        k4q1 = dt * dHdp1(osc.q1+k3q1,osc.L)
        k4q2 = dt * dHdp2(osc.q1+k3q1,osc.L)
        osc.p1 = osc.p1 + 1/6*(k1p1+2*k2p1+2*k3p1+k4p1)
        osc.p2 = osc.p2 + 1/6*(k1p2+2*k2p2+2*k3p2+k4p2)
        osc.q1 = osc.q1 + 1/6*(k1q1+2*k2q1+2*k3q1+k4q1)
        osc.q2 = osc.q2 + 1/6*(k1q2+2*k2q2+2*k3q2+k4q2)
    


# Animation function which integrates a few steps and return a line for the pendulum
def animate(framenr,integrator,pendulum_lines,stepsperframe):
    for it in range(stepsperframe):
        integrator.integrate(osc,obs)

    x1 = math.sin(osc.q1);
    y1 = -math.cos(osc.q1);
    x2 = x1 + math.sin(osc.q2);
    y2 = y1 - math.cos(osc.q2)
    pendulum_lines.set_data([0,x1,x2],[0,y1,y2])
    return pendulum_lines,class Simulation:

    def run(self,tmax=30.,# final time
        stepsperframe=5,# how many integration steps between visualising frames
        outfile='energy1.pdf'
        ):
        numframes = int(tmax / (stepsperframe * integrator.dt))

        plt.clf()
        fig = plt.figure()
        ax = plt.subplot(xlim=(-2.2,2.2),ylim=(-2.2,2.2))
        plt.axhline(y=0)  # draw a default hline at y=1 that spans the xrange
        plt.axvline(x=0)  # draw a default vline at x=1 that spans the yrange
        pendulum_lines,= ax.plot([],[],lw=5)

        oscillator = Oscillator()
        obs = Observables()
        # Call the animator,blit=True means only re-draw parts that have changed
        anim = animation.FuncAnimation(fig,animate,# init_func=init,fargs=[oscillator,stepsperframe],frames=numframes,interval=25,blit=True,repeat=False)
        plt.show()

        # If you experience problems visualizing the animation and/or
        # the following figures comment out the next line
        plt.waitforbuttonpress(30)

        plt.figure()
        plt.xlabel('q1')
        plt.ylabel('p1')
        plt.plot(obs.q1list,obs.p1list)
        plt.tight_layout()  # adapt the plot area tot the text with larger fonts

        plt.figure()
        plt.xlabel('q2')
        plt.ylabel('p2')
        plt.plot(obs.q2list,obs.p2list)
        plt.tight_layout()  # adapt the plot area tot the text with larger fonts

        plt.figure()
        plt.xlabel('q1')
        plt.ylabel('p1')
        plt.plot(obs.poincare_q1,obs.poincare_p1,'ro')
        plt.tight_layout()  # adapt the plot area tot the text with larger fonts

        plt.figure()
        plt.xlabel('time')
        plt.ylabel('energy')
        plt.plot(obs.time,obs.epot,obs.time,obs.ekin,obs.etot)
        plt.legend(('Epot','Ekin','Etot'))
        plt.tight_layout()  # adapt the plot area tot the text with larger fonts
        plt.show()


simulation = Simulation()
simulation.run(integrator=RK4Integrator())

我希望有人可以帮助我解决此问题。预先感谢！

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）