在python中求解SHM的Runge-Kutta 4算法可能的误差来源是什么

Question

我用Python语言写了一个简单的RK4程序来求解SHM方程：

y''[t] = -w^2*y[t]

通过写下：

z'[t] = -w^2*y
y'[t] = z

以下是RK4代码：

def rk4_gen(t, y, z, h):
    k1 = h*f1(y, z, t)
    l1 = h*f2(y, z, t)

    k2 = h*f1(y+k1/2.0, z+l1/2.0, t+h/2.0)
    l2 = h*f2(y+k1/2.0, z+l1/2.0, t+h/2.0)

    k3 = h*f1(y+k2/2.0, z+l2/2.0, t+h/2.0)
    l3 = h*f2(y+k2/2.0, z+l2/2.0, t+h/2.0)

    k4 = h*f1(y+k3, z+l3, t+h)
    l4 = h*f2(y+k3, z+l3, t+h)

    y = y + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0
    z = z + l1/6.0 + l2/3.0 + l3/3.0 + l4/6.0
    t = t + h

    return t, y, z

def f1(y, z, t):
    return z

def f2(y, z, t):
    return -y

结果的曲线图几乎是正弦的，但值超过-1和1的量很小。

Step  0 : t =  0.000, y =   0.00000000, z =   1.00000000
Step  1 : t =  0.100, y =   0.09987500, z =   0.99583542
Step  2 : t =  0.200, y =   0.19891812, z =   0.98171316
Step  3 : t =  0.300, y =   0.29613832, z =   0.95775779
Step  4 : t =  0.400, y =   0.39056108, z =   0.92419231
Step  5 : t =  0.500, y =   0.48123826, z =   0.88133615
Step  6 : t =  0.600, y =   0.56725756, z =   0.82960208
Step  7 : t =  0.700, y =   0.64775167, z =   0.76949228
Step  8 : t =  0.800, y =   0.72190710, z =   0.70159347
Step  9 : t =  0.900, y =   0.78897230, z =   0.62657115
Step  10 : t =  1.000, y =   0.84826536, z =   0.54516314
Step  11 : t =  1.100, y =   0.89918085, z =   0.45817226
Step  12 : t =  1.200, y =   0.94119609, z =   0.36645847
Step  13 : t =  1.300, y =   0.97387644, z =   0.27093037
Step  14 : t =  1.400, y =   0.99687982, z =   0.17253614
Step  15 : t =  1.500, y =   1.00996028, z =   0.07225423
Step  16 : t =  1.600, y =   1.01297061, z =  -0.02891646
Step  17 : t =  1.700, y =   1.00586398, z =  -0.12996648
Step  18 : t =  1.800, y =   0.98869457, z =  -0.22988588
Step  19 : t =  1.900, y =   0.96161722, z =  -0.32767438
Step  20 : t =  2.000, y =   0.92488601, z =  -0.42235127
Step  21 : t =  2.100, y =   0.87885191, z =  -0.51296534
Step  22 : t =  2.200, y =   0.82395944, z =  -0.59860439
Step  23 : t =  2.300, y =   0.76074238, z =  -0.67840440
Step  24 : t =  2.400, y =   0.68981857, z =  -0.75155827
Step  25 : t =  2.500, y =   0.61188388, z =  -0.81732398
Step  26 : t =  2.600, y =   0.52770540, z =  -0.87503206
Step  27 : t =  2.700, y =   0.43811390, z =  -0.92409250
Step  28 : t =  2.800, y =   0.34399560, z =  -0.96400066
Step  29 : t =  2.900, y =   0.24628344, z =  -0.99434256
Step  30 : t =  3.000, y =   0.14594781, z =  -1.01479910
Step  31 : t =  3.100, y =   0.04398694, z =  -1.02514942
Step  32 : t =  3.200, y =  -0.05858305, z =  -1.02527330
Step  33 : t =  3.300, y =  -0.16073825, z =  -1.01515248
Step  34 : t =  3.400, y =  -0.26145719, z =  -0.99487106
Step  35 : t =  3.500, y =  -0.35973108, z =  -0.96461480
Step  36 : t =  3.600, y =  -0.45457385, z =  -0.92466944
Step  37 : t =  3.700, y =  -0.54503210, z =  -0.87541801
Step  38 : t =  3.800, y =  -0.63019464, z =  -0.81733718
Step  39 : t =  3.900, y =  -0.70920170, z =  -0.75099262
Step  40 : t =  4.000, y =  -0.78125355, z =  -0.67703353
Step  41 : t =  4.100, y =  -0.84561868, z =  -0.59618627
Step  42 : t =  4.200, y =  -0.90164114, z =  -0.50924723
Step  43 : t =  4.300, y =  -0.94874724, z =  -0.41707502
Step  44 : t =  4.400, y =  -0.98645148, z =  -0.32058195
Step  45 : t =  4.500, y =  -1.01436144, z =  -0.22072502
Step  46 : t =  4.600, y =  -1.03218196, z =  -0.11849644
Step  47 : t =  4.700, y =  -1.03971818, z =  -0.01491378
Step  48 : t =  4.800, y =  -1.03687770, z =   0.08899018
Step  49 : t =  4.900, y =  -1.02367164, z =   0.19217774
Step  50 : t =  5.000, y =  -1.00021473, z =   0.29361660

为什么y或z的值超出范围[-1,1]？是否有任何正在传播的type-casting错误？？

Answer 1

不，它要么是步长，要么就是浮点数的工作方式。

将步长减半，看看误差是否会减小。

你需要明白，所有像这样的数值算法都是近似值。您不应该期望在正弦波的顶部看到完美的1.0，或者当它穿过x轴时看到完美的0.0。

你不能用二进制精确地表示0.1，就像你不能精确地用十进制表示1/3一样。然后就是IEEE floating point numbers work的方式。你需要理解这些东西。

这里有另一个想法:您对RK4实现是否正确有多大信心？我现在没有时间仔细研究它，但我想知道为什么您不使用像NumPy这样的库，而使用您自己的库？也许他们比我们更了解数值方法，更多的用户意味着更快地发现和修复bug。

你在评论中写了TDSE。我假设这就是Time Dependent Schrodinger Equation，对吗？如果是，您是否在解决方案中正确使用了复数？集成方案是否正确对待它们？