利用Python绘画双摆操作分享
更新时间：2022年04月01日 08:36:11 作者：微小冷
这篇文章主要介绍了利用Python画双摆，绘画双摆的过程主要包括以下步骤，双摆问题、运动过程及公式推导过程，下文详细介绍，需要的小伙伴可以参考一下
1.双摆问题
所谓双摆，就是两个连在一起的摆。
接下来本来是要推公式的，考虑考虑到大家可能会有公式恐惧症，同时又喜欢看图，所以把公式挪到后面。
所以，只需知道角速度的微分方程，就可写出对应的代码，其方程如下:
从而转为代码得到：
# 其中，lam,mu,G_L1,M为全局变量 def derivs(state, t): dydx = np.zeros_like(state) th1,om1,th2,om2 = state dydx[0] = state[1] delta = state[2] - state[0] cDelta, sDelta = np.cos(delta), np.sin(delta) sTh1,_,sTh2,_ = np.sin(state) den1 = M - mu*cDelta**2 dydx[1] = (mu * om1**2 * sDelta * cDelta + mu * G_L1 * sTh2 * cDelta + mu * lam * om2**2 * sDelta - M * G_L1 * sTh1)/ den1 dydx[2] = state[3] den2 = lam * den1 dydx[3] = (- mu * lam * om2**2 * sDelta * cDelta + M * G_L1 * sTh1 * cDelta - M * om1**2 * sDelta - M * G_L1 * sTh2)/ den2 return dydx
接下来根据微分方程的解，便可进行绘图。
# 这段代码用于设置初值，并调用integrate求解微分方程组 import numpy as np import scipy.integrate as integrate G = 9.8 L1,L2 = 1.0, 1.0 G_L1 = G/L1 lam = L2/L1 #杆长度比L2/L1 mu = 1.0 #质量比M2/M1 M = 1+mu # 生成时间 dt = 0.01 t = np.arange(0, 20, dt) th1,th2 = 120.0, -10.0 #初始角度 om1,om2 = 0.0, 0.00 #初始角速度 state = np.radians([th1, om1, th2, om2]) # 微分方程组数值解 y = integrate.odeint(derivs, state, t) # 真实坐标 x1 = L1*sin(y[:, 0]) y1 = -L1*cos(y[:, 0]) x2 = L2*sin(y[:, 2]) + x1 y2 = -L2*cos(y[:, 2]) + y1
至此，就得到了所有位置处的坐标，从而可以观察到双摆的轨迹如图所示
绘图代码为：
import matplotlib.pyplot as plt plt.scatter(x1,y1,marker='.') plt.scatter(x2,y2,marker='.') plt.show()
若将时间设置得长一点，然后在画图的时候更改一下颜色，就会看到双摆的运动区间，可见自然界还是挺有情怀的
其绘图代码为：
plt.plot(x1,y1,marker='.',alpha=0.2, linewidth=0.2)
plt.plot(x2,y2,marker='.',alpha=0.2, linewidth=2, c='r')
plt.axis('off')
plt.show()当然，也可以将其运动轨迹以一种三维的形式绘制出来
ax = plt.gca(projection='3d') ax.plot3D(t,x1,y1,linewidth=1) plt.show()
额……好吧，看来并没有什么情怀。
但是，如果把这两个小球分别当作两个星球，而我们又在一颗星球上，那么所观测到的另一颗星球的运动大致如下，不出意外是个圆，毕竟圆形二者之间的距离是恒定的。
绘图代码为：
ax = plt.gca(projection='3d') ax.plot3D(t,x2-x1,y2-y1,linewidth=0.5) plt.show()
如果更改一下初值，则图形将有如下变化
初值设为：
th1,th2 = 0, 0 #初始角度 om1,om2 = 120.0, 108.00 #初始角速度
2.运动过程
最后，还是传统技能，绘制一下双摆的运动过程如下：
代码为：
import matplotlib.animation as animation
# 下面为绘图过程
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect('equal')
ax.grid()
line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
# 初始化图形
def init():
line.set_data([], [])
time_text.set_text('')
return line, time_text
def animate(i):
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]
line.set_data(thisx, thisy)
time_text.set_text(time_template % (i*dt))
return line, time_text
ani = animation.FuncAnimation(fig, animate, range(1, len(y)),
interval=dt*1000, blit=True, init_func=init)
ani.save("dua_1.gif",writer='imagemagick')
plt.show()3.公式推导过程
双摆的动能和势能分别为:
根据拉格朗日方程
则有：
其中，
展开可得则:
到此这篇关于利用Python画双摆的文章就介绍到这了,更多相关Python画双摆内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！
相关文章
这篇文章主要为大家详细介绍了python实现飞机大战游戏，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下2019-11-11
这篇文章主要介绍了Python数据结构与算法之二叉树结构定义与遍历方法,结合实例形式详细分析了Python实现二叉树结构的定义、遍历方法及相关注意事项,需要的朋友可以参考下2017-12-12
在PyCharm中安装Git时遇到问题,按照视频步骤操作后发现没有git选项,重新检查设置,发现git目录配置错误,重新选择正确的目录后,通过Test确认无误,在使用commit提交时遇到错误,按刷新按钮即可解决2024-11-11
本文给大家介绍三种判断文件或文件夹是否存在的方法，分别使用os模块、Try语句、pathlib模块。感兴趣的朋友一起看看吧2017-07-07
这篇文章主要介绍了python 有效的括号的实现代码示例，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧2019-11-11
这篇文章主要介绍了python初步实现word2vec操作，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧2020-06-06
Ollama提供了Python SDK,使得开发者能够在Python环境中轻松集成和使用本地运行的模型进行自然语言处理任务,具有一定的参考价值,感兴趣的可以了解一下2025-03-03
这篇文章主要介绍了Python如何提取chm数据，具有很好的参考价值，希望对大家有所帮助。如有错误或未考虑完全的地方，望不吝赐教2023-01-01
本文主要介绍了Pycharm连接远程mysql报错的实现,文中通过图文介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧2023-08-08
这篇文章主要介绍了浅析Python中else语句块的使用技巧的相关资料,非常不错，具有参考借鉴价值，需要的朋友可以参考下2016-06-06
最新评论