python - 摆的数值解
问题描述
我在提供我的系统的图形表示时遇到了麻烦,它恰好是一个谐波驱动的钟摆。问题如下图供参考。 问题
我使用的源代码如下所示,使用 Verlet 方案。
#Import needed modules
import numpy as np
import matplotlib.pyplot as plt
#Initialize variables (Initial conditions)
g = 9.8 #Gravitational Acceleration
L = 2.0 #Length of the Pendulum
A0 = 3.0 #Initial amplitude of the driving acceleration
v0 = 0.0 #Initial velocity
theta0 = 90*np.pi/180 #Initial Angle
drivingPeriod = 20.0 #Driving Period
#Setting time array for graph visualization
tau = 0.1 #Time Step
tStop = 10.0 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+tau, tau) #Array of time
theta = np.zeros(len(t))
v = np.zeros(len(t))
#Verlet Method
theta[0] = theta0
v[0] = v0
for i in range(len(t)-1):
accel = -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta[i])
theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]
v[i+1] = v[i] + 0.5*tau*(accel[i] + accel[i+1])
#Plotting and saving the resulting graph
fig, ax1 = plt.subplots(figsize=(7.5,4.5))
ax1.plot(t,theta*(180/np.pi))
ax1.set_xlabel("Time (t)")
ax1.set_ylabel("Theta")
plt.show()
显示了一个示例输出。 输出
钟摆应该刚刚回到它的初始角度。我该如何解决这个问题?请注意,随着时间的推移,我的角度测量值(度)也会增加。我希望它只有 0 度到 360 度的域。
解决方案
请更改数组计算
accel = -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta[i])
theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]
在正确的位置进行正确的元素计算
theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]
accel[i+1] = -((g + (A0*np.sin((2*np.pi*t[i+1]) / drivingPeriod)))/L) * np.sin(theta[i+1])
请注意,您需要accel[0]
在循环外单独计算。
如果您将物理模型的细节分开并在开始时声明,它会使代码更具可读性
def acceleration(t,theta):
return -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta)
这样以后你就可以打电话了
accel[i+1]=acceleration(t[i+1],theta[i+1])
即便如此,通过强制振荡,您的系统是打开的,强制作用可能会将能量泵入钟摆,使其开始旋转。这就是您的图表显示的内容。
与任何辛方法一样,Verlet 方法仅在系统封闭且保守的情况下承诺在某种程度上具有恒定能量,也就是说,在最常见的情况下,没有外部影响并且所有力都是梯度力。
推荐阅读
- python - 我有两个 csv 文件需要使用 for 循环相互比较
- vue.js - 浏览器后退按钮在路由器查询参数更改时无法正常工作
- java - 如何在java中使用缓存管理器获取所有缓存键和值
- javascript - 范围错误 超出最大调用堆栈大小 Angular
- python - 用 csv 文件中的数据替换硬编码代码
- javascript - 我如何跟踪按键并将它们记录到 Javascript 中的控制台
- javascript - 在运行时监控和记录 JavaScript/TypeScript 中的 NaN 值?
- javascript - 将 div 直接移动到另一个带有普通/原生 js 的 div 之后?
- ffmpeg - 如何控制 ffmpeg 中的 mjpeg 压缩质量?
- mysql - 如果在 MySQL 中满足特定条件,如何从一列中减去?