钟摆的数值解

如何解决钟摆的数值解

我无法提供我的系统的图形表示，这恰好是一个谐波驱动的钟摆。问题如下图，供参考。 Problem

我使用的源代码如下所示，使用 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()

显示了示例输出。 Output

钟摆应该回到其初始角度。我该如何解决这个问题？请注意，随着时间的推移，我的角度测量（度）也会增加。我希望它只有 0 度到 360 度的域。

解决方法

请更改数组计算

<span>0</span>
<span>0</span>
<span class="no-bg">:</span>
<span>0</span>
<span>0</span>

在正确的位置进行适当的元素计算

    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)

---

即便如此，在强制振荡的情况下，您的系统是开放的，强制作用可能会将能量泵入钟摆，使其开始旋转。这就是您的图表显示的内容。

与任何辛方法一样，Verlet 方法只承诺在系统是封闭和保守的情况下具有恒定的能量，也就是说，在最常见的情况下，没有外部影响，所有力都是梯度力。