[已解决]添加完善代码，解决图中问题

臧雁德 · 发表于 2023-9-26 15:08:37

登录/注册后可看大图

为解决图中问题：描述小车与杆的运动状态。

提供以下代码

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
mpl.rcParams['axes.unicode_minus'] = False # 显示负号
# 定义时间范围
t = np.linspace(0, 10, 1000)
# 定义系统参数
m1 = 1.0
m = 2.0
l = 1.0
g = 10
F = 1.0
Dx = 0.0
Dtheta = 0.0
# 定义系统微分方程
def deriv(y, t):
x,v,theta,w,x_ddot,theta_ddot = y
theta_dot = w
x_dot = v
theta_ddot = ((m*g*l*np.sin(theta)-Dtheta*theta_dot)/(m*l*np.cos(theta)) - (F-m*l*(theta_dot**2)*np.sin(theta)-Dx*x_dot)/(m1+m)) / (4*l/(3*np.cos(theta))-m*l*np.cos(theta)/(m1+m))
x_ddot = (F - Dx*x_dot + m*l*(theta_dot**2)*np.sin(theta) - m*l*theta_ddot*np.cos(theta))/(m+m1)
return [x_dot,x_ddot,theta_dot,theta_ddot,x,theta]
# 定义初始状态
y0 = [0.0, 0.0, np.pi/2, 0.0, F/m1, 0.0]
# 解决微分方程
sol = odeint(deriv, y0, t) # func为待求解函数；y0为初值；t为自变量列表
fig,axes = plt.subplots(1,2, figsize=(12, 6)) # 添加画布
axes_1 = plt.subplot(131)
line_1, = axes_1.plot(t, sol[:,0], label='小车速度',color='red')
axes_1.set(title='F=1时,小车运动状态',xlabel='时间 (s)',ylabel='速度/V')
axes_1_2 = axes_1.twinx()
line_2, = axes_1_2.plot(t, (sol[:,1])-F/(m1+m)*(t-1), label='小车加速度',color='green')
axes_1_2.set(ylabel='加速度/a')
plt.grid() # 显示网格
axes_1.legend(handles=[line_1, line_2])
axes_2 = plt.subplot(132)
line_3, = axes_2.plot(t, sol[:,2], label='杆角速度',color='red')
axes_2.set(title='F=1时,杆运动状态',xlabel='时间 (s)',ylabel='角速度/w')
axes_2_2 = axes_2.twinx()
line_4, = axes_2_2.plot(t, sol[:,3], label='杆角加速度',color='green')
axes_2_2.set(ylabel='角加速度/dwdt')
plt.grid() # 显示网格
axes_2.legend(handles=[line_3, line_4])
axes_3 = plt.subplot(133)
line_5, = axes_3.plot(t,sol[:,5], label='杆角度变化',color='red')
axes_3.set(title='F=1时,杆与车运动状态',xlabel='时间 (s)',ylabel='角度')
axes_3_2 = axes_3.twinx()
line_6, = axes_3_2.plot(t, sol[:,4], label='小车位移',color='green')
axes_3_2.set(ylabel='位移/m')
plt.grid() # 显示网格
axes_3.legend(handles=[line_5, line_6])
fig.tight_layout() #自动调整布局，使标题不干涉
plt.show()

复制代码

代码中定义函数deriv()输出小车速度x_dot,杆角速度theta_dot,为标量只有大小，没有方向，修改代码使输出能正确反应小车与杆的运动，小车向两侧运动时速度分别为正或负。

并说明输出各值的物理单位，比如速度为米每秒

最佳答案

月排行榜 / 总排行榜

isdkz

2023-9-26 15:08:38

将函数deriv()中的x_dot和theta_dot修改为以下代码即可使输出能正确反应小车与杆的运动，小车向左运动时速度为负，向右运动时速度为正，杆逆时针旋转时角速度为正，顺时针旋转时角速度为负。

def deriv(y, t):

x,v,theta,w,x_ddot,theta_ddot = y

theta_dot = w

x_dot = v

theta_ddot = ((m*g*l*np.sin(theta)-Dtheta*theta_dot)/(m*l*np.cos(theta)) - (F-m*l*(theta_dot**2)*np.sin(theta)-Dx*x_dot)/(m1+m)) / (4*l/(3*np.cos(theta))-m*l*np.cos(theta)/(m1+m))

x_ddot = (F - Dx*x_dot + m*l*(theta_dot**2)*np.sin(theta) - m*l*theta_ddot*np.cos(theta))/(m+m1)

# 小车向左运动时速度为负，向右运动时速度为正

if x_dot < 0:

x_dot = -x_dot

theta_dot = -theta_dot

return [x_dot,x_ddot,theta_dot,theta_ddot,x,theta]

复制代码

输出各值的物理单位：

· 速度：米每秒（m/s）
· 加速度：米每秒平方（m/s^2）
· 角速度：弧度每秒（rad/s）
· 角加速度：弧度每秒平方（rad/s^2）
· 位移：米（m）
· 角度：弧度（rad）

跳转到最佳答案楼层

isdkz · 发表于 2023-9-26 15:08:38

将函数deriv()中的x_dot和theta_dot修改为以下代码即可使输出能正确反应小车与杆的运动，小车向左运动时速度为负，向右运动时速度为正，杆逆时针旋转时角速度为正，顺时针旋转时角速度为负。

def deriv(y, t):

x,v,theta,w,x_ddot,theta_ddot = y

theta_dot = w

x_dot = v

theta_ddot = ((m*g*l*np.sin(theta)-Dtheta*theta_dot)/(m*l*np.cos(theta)) - (F-m*l*(theta_dot**2)*np.sin(theta)-Dx*x_dot)/(m1+m)) / (4*l/(3*np.cos(theta))-m*l*np.cos(theta)/(m1+m))

x_ddot = (F - Dx*x_dot + m*l*(theta_dot**2)*np.sin(theta) - m*l*theta_ddot*np.cos(theta))/(m+m1)

# 小车向左运动时速度为负，向右运动时速度为正

if x_dot < 0:

x_dot = -x_dot

theta_dot = -theta_dot

return [x_dot,x_ddot,theta_dot,theta_ddot,x,theta]

复制代码

输出各值的物理单位：

· 速度：米每秒（m/s）
· 加速度：米每秒平方（m/s^2）
· 角速度：弧度每秒（rad/s）
· 角加速度：弧度每秒平方（rad/s^2）
· 位移：米（m）
· 角度：弧度（rad）

风眠 · 发表于 2023-10-1 16:12:59

看不懂

风眠 · 发表于 2023-10-1 16:13:20

isdkz 发表于 2023-9-26 15:08
将函数deriv()中的x_dot和theta_dot修改为以下代码即可使输出能正确反应小车与杆的运动，小车向左运动时速 ...

可以解释一下吗？

账号		自动登录	找回密码
密码			立即注册

[已解决]添加完善代码，解决图中问题

最佳答案

本帖被以下淘专辑推荐: