I have this piece of code (I put just a part of it) in 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])
So what I am doing is to solve the simple harmonic oscillator using a symplectic Runge-Kutta. (The formulas I use are from a book and I checked I copied them correctly). So in iterate_Z I pick some initial values for Z and I do some iterations, in order to get a much better value, from which to go to the next step of Runge-Kutta. When I run the program (and add a print inside the function), it gives output for a few steps and after that I get this error:
File "mtrand.pyx", line 2093, in mtrand.RandomState.f (numpy/random/mtrand/mtrand.c:23649)ValueError: dfnum <= 0
pointing to this line:
sum = sum + a[i][j]*f(y[-1]+Z[j],t[-1])
Can someone tell me what this error means and how can I fix it? Thank you!
