作为数值分析技术的新手，尤其是 RK2，我认为最好的方法是使用 python 来解决众所周知的使用 RK2 技术的质量弹簧振荡器。

我的问题是步长似乎会影响我的解决方案的周期，我知道为什么，如果有人能帮助我更好地理解我做错了什么，我将不胜感激。欢迎评论代码！谢谢你。

问题如下（这不是教科书或家庭作业问题！）： 质量块连接到具有线性弹簧常数的弹簧。当时$t=0$质量偏离平衡，$x=0$, 到位置$x(t=0) = 1 cm$. 质量被释放并在平衡附近振荡。使用 RK2 的方法绘制质量位置随时间的演变。该问题所需的一阶微分方程为：

$\dot{x} = -\omega x(t)$

$\frac{dx}{dt} = \dot{x}$

我使用的 RK2 功能如下：

$v(t_{n+1}) = v(t_n) -\omega x(t_n)\delta t$

$x(t_{n+1}) = x(t_n) + v(t_n)\delta t$

$k2_{x} = \delta t (v(t_n) - \omega x(t_n)\frac{\delta t}{2})$

$k2_{v} = \delta t (\omega x(t_n) + v(t_n) \frac{\delta t}{2})$

最后两个函数是我认为用于 RK2 方法的中点值。

连同实现的代码：

import numpy as np
import matplotlib.pyplot as plt

print "RK2 Method for Oscillating Spring"

#Spring Constant
k = 1.0
#Mass
m = 0.02
#Angular frequency squared
w = k/m

period = 2*np.pi*(np.sqrt(1/w))

#Initial Conditions
x_start = 1.0
v_start = 0.0

#Set loop variables
x_position = x_start
velocity = v_start

#Step Size
dt = 0.001

#Arrays to store variable values
x_plot = []
v_plot = []
time = []

#Plot the numerically determined position over
# the time interval 1 to 10 seconds in steps of 0.1
for t in np.arange(1,10,0.1):

    #The two equations below are used in the Euler Method approach
    #Their solution's period, which can be plotted, varies depending on the step size
    # and I know it shouldn't
    #x_position = x_position + dt*velocity
    #velocity = velocity - dt*w*x_position

    #The equations below are used for RK2 method. I don't know why I cant get it to work.
velocity = velocity + dt*(-w*x_position)
    k2_x = dt*(velocity + (-w*x_position*dt/2.0))
    x_position = x_position + k2_x

    k2_v = -dt*w*(x_position + velocity*dt/2.0)
    velocity = velocity + k2_v


    x_plot.append(x_position)
    v_plot.append(velocity)
    time.append(t)


plt.subplot(3,1,1)
plt.plot(time,x_plot)
plt.subplot(3,1,2)
plt.plot(time,v_plot)
plt.subplot(3,1,3)
plt.plot(x_plot,v_plot) # illustrate divergence over number of steps
plt.show()

我的问题是我不明白为什么改变步长 dt 会改变解的周期（正弦函数）。