Ошибка конвергенции для Рунге-Кутты в Python
У меня есть этот кусок кода (я положил только часть этого) в Python:
from __future__ import division
import numpy as np
from pylab import *
from numpy import linalg as LA
a = [[5/36,2/9-sqrt(15)/15,5/36-sqrt(15)/30],[5/36+sqrt(15)/24,2/9,5/36-sqrt(15)/24],[5/36+sqrt(15)/30,2/9+sqrt(15)/15,5/36]]
treshold = 10**(-12)
h = 0.01
def sho(x, t=None):
k = 1.
t = float(t)
x = array(x)
return array([ x[1],
- x[0]])
def iterate_Z(func,Z,y,h,t):
crit = 1
Z = transpose(Z)
while crit > treshold:
T = []
for i in range(len(Z)):
sum = 0
for j in range(len(Z)):
sum = sum + a[i][j]*f(y[-1]+Z[j],t[-1])
T = T + [h*sum]
crit = norm(Z-T)
Z = np.array(T)
return(Z)
iterate_Z(sho,[[0,0,0],[0,0,0]],[[1,2],[2,3]],0.01,[0,1])
Итак, что я делаю, так это решаю простой гармонический осциллятор с помощью симплектической Рунге-Кутты. (Формулы, которые я использую, взяты из книги, и я проверил, правильно ли я их скопировал). Поэтому в iterate_Z я выбираю некоторые начальные значения для Z и делаю несколько итераций, чтобы получить гораздо лучшее значение, из которого можно перейти к следующему шагу Рунге-Кутта. Когда я запускаю программу (и добавляю печать внутри функции), она выдает вывод за несколько шагов, и после этого я получаю эту ошибку:
File "mtrand.pyx", line 2093, in mtrand.RandomState.f (numpy/random/mtrand/mtrand.c:23649)ValueError: dfnum <= 0
указывая на эту строку:
sum = sum + a[i][j]*f(y[-1]+Z[j],t[-1])
Может кто-нибудь сказать мне, что означает эта ошибка и как я могу это исправить? Спасибо!