提问人：Elisa 提问时间：11/7/2023 最后编辑：phoElisa 更新时间：11/8/2023 访问量：53
在 Python 上解决 SHM 问题
Solving SHM on Python
问：
我需要使用欧拉方法迭代求解一个简单的谐波振荡器系统。这是我到目前为止的代码。
#importing necessary libraries
import numpy as np
import matplotlib.pyplot as plt
#defining constants
x_0 = 1e-10 #initial position
v_0 = 0 #initial velocity
m_sodium = 3.82e-26 #mass of sodium atom
k = 12.2 #force constant
#function for potential energy
def potential(x):
Vx = k * (x ** 2) / 2
return Vx
#equation angular frequency
omega = np.sqrt(k / m_sodium)
#creating an array for time
t = np.linspace(0, 10, 1000)
#initialising x and v arrays (to eventually hold the results)
x = np.zeros(len(t))
v = np.zeros(len(t))
#initial conditions
x[0] = x_0
v[0] = v_0
#defining time-step
dt = t[1] - t[0]
#Euler method to solve ODE for all the iterations
for i in range(1, len(t)):
x[i] = v[i] * dt
v[i] = omega**2 * x[i] * dt
#plotting a graph x vs t
plt.plot(t, x)
plt.xlabel('Time t')
plt.ylabel('Position x')
plt.grid(True)
plt.title('Simple Harmonic Oscillator System')
我本来是要获得正弦波的，但是我得到的却是我附加的图像。谁能指出我做错了什么？
答：
2赞
lastchance
11/8/2023
#1
你要解的方程是
DX/DT=V
DV/DT=-欧米茄2 x
如果将这些组合在一起，您将得到标准的 SHM 方程
d 2 x/dt 2 = -欧米茄2x
具有自然循环频率欧米茄和周期 2.pi / 欧米茄
dt=0.01 的时间步长对于您根据该力常数计算出的 omega 来说太大了。相反，您可以将跟踪的整个长度更改为 10 个周期，并使用
t = np.linspace(0, 10 * 2 * np.pi / omega, 1000)
你的更新步骤是错误的（反正一阶欧拉不是很好）。
#importing necessary libraries
import numpy as np
import matplotlib.pyplot as plt
#defining constants
x_0 = 1 #initial position
v_0 = 0 #initial velocity
m_sodium = 3.82e-26 #mass of sodium atom
k = 12.2 #force constant
#equation angular frequency
omega = np.sqrt(k / m_sodium)
#creating an array for time
t = np.linspace(0, 10 * 2 * np.pi / omega, 1000) # 10 oscillation periods
#initialising x and v arrays (to eventually hold the results)
x = np.zeros(len(t))
v = np.zeros(len(t))
#initial conditions
x[0] = x_0
v[0] = v_0
#defining time-step
dt = t[1] - t[0]
#Euler method to solve ODE for all the iterations (THERE ARE BETTER METHODS)
for i in range(1, len(t)):
x[i] = x[i-1] + v[i-1] * dt
v[i] = v[i-1] - omega**2 * x[i] * dt
#plotting a graph x vs t
plt.plot(t, x)
plt.xlabel('Time t')
plt.ylabel('Position x')
plt.grid(True)
plt.title('Simple Harmonic Oscillator System')
plt.show()
输出：
评论
x[1] = v[1] * dtv[1]x[1]x[i] += v[i-1]*dt