Search code examples
pythonnumpyindexoutofboundsexceptionodeintindex-error

I can't solve issue "axis -1 is out of bounds for array of dimension 0"

I'm trying to model the motion of a spring pendulum and this is my code:

import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi, array
import matplotlib.pyplot as plt


#Specify initial conditions
init = array([pi/18, 0]) # initial values

def deriv(z, t):
    x,y=z
    dy=np.diff(y,1)
    dy2=np.diff(y,2)
    dx=np.diff(x,1)
    dx2=np.diff(x,2)
    dt=np.diff(t,1)
    dt2=np.diff(t,1)
    dx2dt2=(4+x)*(dydt)^2-5*x+9.81*cos(y)
    dy2dt2=(-9.81*sin(y)-2*(dxdt)*(dydt))/(l+x)
    return np.array([dx2dt2,dy2dt2])

time = np.linspace(0.0,10.0,1000)
y = odeint(deriv,init,time)

plt.xlabel("time")
plt.ylabel("y")
plt.plot(time, y)
plt.show()

I keep getting the error 

Traceback (most recent call last):
  File "/Users/cnoxon/Desktop/GRRR.py", line 24, in <module>
    y = odeint(deriv,init,time)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/scipy/integrate/odepack.py", line 233, in odeint
    int(bool(tfirst)))
  File "/Users/cnoxon/Desktop/GRRR.py", line 13, in deriv
    dy=np.diff(y,1)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/numpy/lib/function_base.py", line 1163, in diff
    axis = normalize_axis_index(axis, nd)
numpy.core._internal.AxisError: axis -1 is out of bounds for array of dimension 0

I'm a complete beginner to Python so I won't understand most of the terminology, so please bear with me. How do I fix this problem? I'm trying to plot the solutions to the two equations 

dx2dt2=(4+x)*(dydt)^2-5*x+9.81*cos(y)
dy2dt2=(-9.81*sin(y)-2*(dxdt)*(dydt))/(l+x)

but I'm having a lot of trouble. Would someone please explain to me how I should rewrite my code to resolve this?

Thank you!

Solution

  • The problem happens because x and y are integers, not arrays, so you can't do np.diff(y,1).

    But your problem is deeper. Each entry of the y array must fully describe your system, this means that every value needed to compute dx2dt2 and dy2dt2 must be in this vector. So y has to be a list of [x, y, dxdt, dydt]. (Adapt init to correspond to this)

    Then, your deriv function just has to give the derivative of such a vector, which is: [dxdt, dydt, dx2dt2, dy2dt2]. Your deriv function becomes very simple!

    def deriv(z, t):
    x, y, dxdt, dydt = z

    dx2dt2=(4+x)*(dydt)^2-5*x+9.81*cos(y)
    dy2dt2=(-9.81*sin(y)-2*(dxdt)*(dydt))/(l+x)

    return np.array([dxdt, dydt, dx2dt2, dy2dt2])

    And you have two other little errors: use ** instead of ^ in python, and I think you changed a 1 into a l...