带弹簧双摆运动Python模拟代码NameError报错问题求解

问题原因 #

• 变量作用域错误：你在三个动力学函数内部引用的y、z、x都不属于对应函数的入参，也不是函数内定义的局部变量，Python无法识别所以抛出NameError。另外你在代码末尾定义的全局变量x/y/z是调用solve_ivp执行完成后才赋值的，动力学函数运行时这些变量还不存在。
• 求解逻辑错误：你的三个运动方程是强耦合的，theta1、theta2、r的导数互相依赖，不能分开调用三次solve_ivp独立求解，必须把所有状态变量合并为一个状态向量同步求解。

修复后的可运行代码 #

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

g = 9.806 # 重力加速度
l0 = 1 # 弹簧原长
k = 2 # 弹簧劲度系数
OA = 2 # OA段长度
m = 1 # 质点质量

def pendulumDynamics(t, s):
    # 状态向量s定义：[theta1, theta1_dot, theta2, theta2_dot, r, r_dot]
    x1, x2, y1, y2, z1, z2 = s
    
    Fs = -k*(z1 - l0)
    # 各状态导数计算
    x1dot = x2
    x2dot = (Fs*sin(y1 - x1) - m*g*sin(x1))/(m*OA)
    
    y1dot = y2
    y2dot = (-g*sin(y1) - (Fs*cos(y1-x1)*sin(x1))/m + g*cos(y1-x1)*sin(x1) - x2*z1*sin(x1))/z1
    
    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 [x1dot, x2dot, y1dot, y2dot, z1dot, z2dot]

# 初始条件定义
d2r = pi/180
# 初始状态：theta1=30度、theta1_dot=0，theta2=60度、theta2_dot=0，r=1、r_dot=0
s0 = [30*d2r, 0, 60*d2r, 0, 1, 0]
t0 = 0
tf = 10

# 单次求解所有耦合ODE
sol = solve_ivp(pendulumDynamics, [t0, tf], s0, dense_output=True)

t = linspace(t0, tf, 200)
s = sol.sol(t)
# 提取各变量：
theta1 = s[0]
theta1_dot = s[1]
theta2 = s[2]
theta2_dot = s[3]
r = s[4]
r_dot = s[5]

后续操作 #

提取出的6个变量就是你需要的所有状态量，直接用matplotlib对应的绘图接口就可以绘制各变量的变化曲线。

内容的提问来源于stack exchange，提问作者Omor Almamun