在 Python 上解决 SHM 问题

Question

我需要使用欧拉方法迭代求解一个简单的谐振子系统。这是我到目前为止的代码。

#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')

我本打算获得正弦波，但是我获得的是我所附加的图像。有人可以指出我做错了什么吗？

Answer 1

你要解的方程是

dx/dt=v
dv/dt=-omega²x

如果将这些结合起来，您将得到标准 SHM 方程

d²x/dt² = -omega² x

具有自然圆频率 omega 和周期 2.pi / omega

您的时间步长 dt=0.01 对于您根据该力常数计算的欧米伽来说太大了。相反，您可以使用

将迹线的整个长度更改为 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()

输出：