Каков возможный источник ошибки Рунге-Кутта 4 алгоритма в питоне для решения SHM

Question

Я написал простую программу RK4 в Python для решения уравнения 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., 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 в верхней части синусоидальной волны или идеальный 0.0, когда он пересекает ось x.

Вы не можете представить 0,1 точно в двоичном виде больше, чем вы можете представлять 1/3 точно в десятичном формате. И тогда есть путь . Вам нужно понять это.

Вот еще одна мысль: Насколько вы уверены, что ваша реализация RK4 верна? У меня нет времени, чтобы прорваться через него прямо сейчас, но мне было бы интересно, почему вы не используете библиотеку, например, NumPy вместо своей? Возможно, они понимают численные методы лучше, чем мы, и больше пользователей означает поиск и исправление ошибок быстрее.

Вы написали TDSE в своем комментарии. Правильно ли я полагаю, что это Time Dependent Schrodinger Equation? Если да, правильно ли вы используете сложные номера в своем решении? Правильно ли работает схема интеграции?