Ошибка конвергенции для Рунге-Кутты в 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])

Может кто-нибудь сказать мне, что означает эта ошибка и как я могу это исправить? Спасибо!

0 ответов

Другие вопросы по тегам