使用Euler法求解简谐振动SHM异常：未得正弦曲线问题排查

尝试使用Euler迭代法求解简谐振动（Simple Harmonic Oscillator, SHM）系统，编写了如下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')

错误分析 #

你的代码存在两个核心问题，直接导致结果异常：

1. 迭代逻辑完全错误
   Euler法的核心是用前一个时间步的状态计算当前步状态，但你在计算x[i]时，使用了尚未赋值的v[i]（初始为0），且完全没有利用上一步的x[i-1]和v[i-1]，完全违背迭代逻辑。
2. 加速度符号错误
   简谐振动的运动方程是a = -ω²x（回复力与位移方向相反，由F=-kx=ma推导而来），你的代码缺少负号，导致加速度方向错误，运动规律完全偏离预期。

修正后的代码 #

#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] = x[i-1] + v[i-1] * dt  # 基于前一步速度更新当前位置
    v[i] = v[i-1] - omega**2 * x[i-1] * 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()

补充说明 #

修正后会得到近似的正弦曲线，但显式Euler法本身存在能量不守恒的问题，模拟时间越长，振幅会逐渐增大。如果需要更稳定的长期模拟，建议使用Verlet积分法或半隐式Euler法，这类方法在保守系统中的表现更优。

内容的提问来源于stack exchange，提问作者Elisa