0

我想用 python 中的弹簧绘制双摆的运动。我需要绘制 theta1、theta2、r 及其一阶导数。我找到了运动方程，它们是二阶 ODE，因此我将它们转换为一阶 ODE，其中 x1=theta1, x2=theta1-dot, y1=theta2, y2=theta2-dot, z1=r, z2=r 点。这是双摆问题的图片：在此处输入图片描述

这是我的代码：

from scipy.integrate import solve_ivp
from numpy import pi, sin, cos, linspace

g = 9.806 #Gravitational acceleration
l0 = 1 #Natural length of spring is 1
k = 2 #K value for spring is 2
OA = 2 #Length OA is 2
m = 1 #Mass of the particles is 1

def pendulumDynamics1(t, x): #Function to solve for theta-1 double-dot
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]

    Fs = -k*(z1-l0)
    T = m*(x2**2)*OA + m*g*cos(x1) + Fs*cos(y1-x1)

    x1dot = x2
    x2dot = (Fs*sin(y1-x1) - m*g*sin(x1))/(m*OA) # angles are in radians
    return [x1dot,x2dot]

def pendulumDynamics2(t, y): #Function to solve for theta-2 double-dot
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]
    Fs = -k*(z1-l0)
    
    y1dot = y2
    y2dot = (-g*sin(y1) - (Fs*cos(y1-x1)*sin(x1))/m + g*cos(y1-x1)*sin(x1) - x2*z1*sin(x1))/z1
    return [y1dot,y2dot]

def pendulumDynamics3(t, z): #Function to solve for r double-dot (The length AB which is the spring)
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]
    Fs = -k*(z1-l0)
    
    z1dot = z2
    z2dot = g*cos(y1) - Fs/m + (y2**2)*z1 + x2*OA*cos(y1-x1) - (Fs*(sin(y1-x1))**2)/m + g*sin(x1)*sin(y1-x1)
    return [z1dot,z2dot]

# Define initial conditions, etc
d2r = pi/180
x0 = [30*d2r, 0] # start from 30 deg, with zero velocity
y0 = [60*d2r, 0] # start from 60 deg, with zero velocity
z0 = [1, 0] #Start from r=1
t0 = 0
tf = 10

#Integrate dynamics, initial value problem
sol1 = solve_ivp(pendulumDynamics1,[t0,tf],x0,dense_output=True) # Save as a continuous solution
sol2 = solve_ivp(pendulumDynamics2,[t0,tf],y0,dense_output=True) # Save as a continuous solution
sol3 = solve_ivp(pendulumDynamics3,[t0,tf],z0,dense_output=True) # Save as a continuous solution

t = linspace(t0,tf,200) # determine solution at these times
dt = t[1]-t[0]
x = sol1.sol(t)
y = sol2.sol(t)
z = sol3.sol(t)

我的代码中有 3 个函数，每个函数用于求解 x、y 和 z。然后我使用solve_ivp 函数求解x、y 和z。代码中的错误是：

`文件“C:\Users\omora\OneDrive\Dokument\AERO 211\project.py”，第 13 行，在 pendulumDynamics1 y1 = y[0]

NameError：名称'y'未定义`

我不明白为什么说 y 没有定义，因为我在我的函数中定义了它。

4

2 回答 2

1

你的系统是封闭的，没有摩擦，因此可以被拉格朗日或哈密顿形式主义捕获。您有 3 个位置变量，因此是 6 维动态状态，由速度或冲量补充。

让,q_k它们的时间导数和冲量变量, 那么动力学可以通过theta_1, theta_2, rDq_kp_kq_k

def DoublePendulumSpring(u,t,params):
    m_1, l_1, m_2, l_2, k, g = params
    q_1,q_2,q_3 = u[:3]
    p = u[3:]

    A = [[l_1**2*(m_1 + m_2), l_1*m_2*q_3*cos(q_1 - q_2), -l_1*m_2*sin(q_1 - q_2)], 
         [l_1*m_2*q_3*cos(q_1 - q_2), m_2*q_3**2, 0], 
         [-l_1*m_2*sin(q_1 - q_2), 0, m_2]]

    Dq = np.linalg.solve(A,p)

    Dq_1,Dq_2,Dq_3 = Dq
    T1 = Dq_2*q_3*sin(q_1 - q_2) + Dq_3*cos(q_1 - q_2)
    T3 = Dq_1*l_1*cos(q_1 - q_2) + Dq_2*q_3

    Dp = [-l_1*(m_2*Dq_1*T1 + g*(m_1+m_2)*sin(q_1)), 
          l_1*m_2*Dq_1*T1 - g*m_2*q_3*sin(q_2), 
          m_2*Dq_2*T3 + g*m_2*cos(q_2) + k*(l_2 - q_3) ]
    return [*Dq, *Dp]

有关推导，请参见欧拉-拉格朗日方程及其与汉密尔顿方程的联系。你可能会被问到这样的推导。

在适当地定义参数元组和初始条件之后，可以将其提供给odeint并产生一个解决方案，然后可以绘制、动画或以其他方式检查。较低的 bob 追踪一条类似于下面的路径，不是周期性的，也不是非常确定的。（上bob的支点和弧线也插入了，但没那么有趣。）

在此处输入图像描述

于 2021-12-05T17:34:03.387 回答
0 
def pendulumDynamics1(t, x):
    x1 = x[0]
    x2 = x[1]
    y1 = y[0]
    y2 = y[1]
    z1 = z[0]
    z2 = z[1]

x您只能作为参数传递。函数内部的代码不知道y和z引用什么。

您将需要更改函数调用以包含这些变量。

def pendulumDynamics1(t, x, y, z):
于 2021-12-05T05:52:39.967 回答