What is the possible source of error in Runge-Kutta 4 algorithm in python for solving SHM

Question

I wrote a simple RK4 program in python to solve SHM equation:

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

by writing:

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

Following is the RK4 code:

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

The plot of the result is coming almost sinusoidal, but the values exceeding beyond -1 and 1 by small amount.

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

Why are the value of y or z exceeding beyond range [-1,1]? Is there any type-casting error that is propagating??

Answer 1

No, it's either step size or just the way floating point numbers work.

Cut the step size in half and see if the error goes down.

You need to understand that all numerical algorithms like this are approximations. You should not expect to see a perfect 1.0 at the top of the sine wave or a perfect 0.0 when it crosses the x-axis.

You can't represent 0.1 exactly in binary any more than you can represent 1/3 exactly in decimal. And then there's the way IEEE floating point numbers work. You need to understand these things.

Here's another thought: How confident are you that your RK4 implementation is correct? I don't have the time to pore through it right now, but I'd wonder why you don't use a library like NumPy instead of your own? Perhaps they understand numerical methods better than we do, and more users means finding and fixing bugs faster.

You wrote TDSE in your comment. Am I correct in assuming that this is the Time Dependent Schrodinger Equation? If yes, are you using complex numbers correctly in your solution? Is the integration scheme treating them properly?