python中的双摆动画

问题描述

我使用以下代码对微分方程组进行数值求解，然后对其进行动画处理：

for i in range(steps-1):
    Theta1 = t1[i]
    Theta2 = t2[i]
    dTheta1 = w1[i]
    dTheta2 = w2[i]
    a1 = (g*(np.sin(Theta2)*np.cos(Theta1-Theta2)-mu*np.sin(Theta1))-(l2*dTheta2*dTheta2+l1*dTheta1*dTheta1*np.cos(Theta1-Theta2))*np.sin(Theta1-Theta2))/(l1*(mu-np.cos(Theta1-Theta2)*np.cos(Theta1-Theta2)))
    a2 = (mu*g*(np.sin(Theta1)*np.cos(Theta1-Theta2)-np.sin(Theta2))+(mu*l1*dTheta1*dTheta1+l2*dTheta2*dTheta2*np.cos(Theta1-Theta2))*np.sin(Theta1-Theta2))/(l2*(mu-np.cos(Theta1-Theta2)*np.cos(Theta1-Theta2)))
    w1[i+1] = w1[i] + dt*a1
    w2[i+1] = w2[i] + dt*a2
    t1[i+1] = t1[i] + dt*w1[i]

    t2[i+1] = t2[i] + dt*w2[i] 

这给了我错误，例如 运行时警告：longdouble_scalars 中遇到无效值 和 RuntimeWarning: longdouble_scalars 中遇到溢出。 起初，我认为这可能是由于精度错误造成的，所以我尝试使用 longdouble。但这根本没有帮助。当我运行动画时，我得到了几帧非物理系统，然后它关闭了。我已经从多个来源检查了我的方程式，最终我也使用了他们的方程式，但问题仍然存在。我做错了什么？

编辑：好的，所以简单地减小步长似乎不再给我错误，但现在动画是非物理的。随着时间的推移，钟摆似乎会加快速度，并且还会进行完整的旋转，即使它从极端静止状态开始也是如此。我正在使用四阶 RK 方法，但仍然遇到此问题。有任何想法吗？我在下面提供了我的新代码

def diff_eq_a1(theta1,theta2,w1,w2):
    a11 = mu
    a12 = np.cos(theta1-theta2)
    a21 = np.cos(theta1-theta2)
    a22 = 1
    
    b1 = (mu*g*l1*np.sin(theta1)) + (l2*np.sin(theta1-theta2)*w2*w2)
    b2 = (g*l2*np.sin(theta2))-(l1*np.sin(theta1-theta2)*w1*w1)
    
    return -(1/l1)*(b1*a22-b2*a12)/(a11*a22-a21*a12)
    

def diff_eq_a2(theta1,w2):
    a11 = mu
    a12 = np.cos(theta1-theta2)
    a21 = np.cos(theta1-theta2)
    a22 = 1
    
    b1 = (mu*g*l1*np.sin(theta1)) + (l2*np.sin(theta1-theta2)*w2*w2)
    b2 = (g*l2*np.sin(theta2))-(l1*np.sin(theta1-theta2)*w1*w1)
    
    return -(1/l2)*(b2*a11-b1*a21)/(a11*a22-a21*a12)


for i in range(steps-1):
    a1_k1 = dt*diff_eq_a1(t1[i],t2[i],w1[i],w2[i])
    a1_k2 = dt*diff_eq_a1((t1[i]+0.5*a1_k1),w2[i])
    a1_k3 = dt*diff_eq_a1((t1[i]+0.5*a1_k2),w2[i])
    a1_k4 = dt*diff_eq_a1((t1[i]+a1_k3),w2[i])
    
    w1[i+1] = w1[i] + (a1_k1/6)+(a1_k2/3)+(a1_k3/3)+(a1_k4/6)
    
    a2_k1 = dt*diff_eq_a1(t1[i],w2[i])
    a2_k2 = dt*diff_eq_a1(t1[i],(t2[i]+0.5*a2_k1),w2[i])
    a2_k3 = dt*diff_eq_a1(t1[i],(t2[i]+0.5*a2_k2),w2[i])
    a2_k4 = dt*diff_eq_a1(t1[i],(t2[i]+a2_k3),w2[i])
    
    w2[i+1] = w2[i] + (a2_k1/6)+(a2_k2/3)+(a2_k3/3)+(a2_k4/6)
    
    t1[i+1] = t1[i] + dt*w1[i]
    t2[i+1] = t2[i] + dt*w2[i]

x1 = l1*np.sin(t1)
x2 = x1 + l2*np.sin(t2)
y1 = - l1*np.cos(t1)
y2 = y1 - l2*np.cos(t2)

我知道它现在不是很整洁，但我真的只想先让它工作

解决方法

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

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

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