Search code examples
pythonnumpyjupyter-notebookrunge-kutta

Index out of bounds in Runge-Kutta algorithm for parabolic shot


Hi I am trying to get this program to run. The z structure is assumed to be size 4, however it gives the error. I also notice that if I print z in 1 it gives me a value, so I do not know what is the flaw. It is a program that calculates, using runge kutta, a physical friction parabolic shot system, however the following error occurs: Image of the error

The code is as follows:

import numpy as np
import matplotlib.pyplot as plt


def RK4(f_user,U0,t,*params):
    n=t.size
    neq=U0.size
    u=np.zeros((n,neq))
    u[0]=U0
    dt=t[2]-t[0]
    
    for i in range(n-1):
        k1=f_user(u[i],t[i],*params)
        k2=f_user(u[i]+k1*dt/2.,t[i]+dt/2.,*params)
        k3=f_user(u[i]+k2*dt/2.,t[i]+dt/2.,*params)
        k4=f_user(u[i]+k3*dt,t[i+1],*params)
        u[i+1]=u[i]+dt*(k1+2*k2+2*k3+k4)/6.
    return u

g = 9.81
# m = 1  # kg
theta_0 = 45 * np.pi / 180  # radianes
v_0 = 10  # m/s
tInc = 0.01
tStop = (1.+ tInc)*20
t = np.arange(0., tStop, tInc)
x_0=0
y_0=0
vx_0=v_0 * np.cos(theta_0)
vy_0=v_0 * np.sin(theta_0)
z0 = np.array([x_0,vx_0,y_0,vy_0])

def f(z,t,N,k,x_0,vx_0,y_0,vy_0):
        v=np.sqrt(z[1] ** 2 + z[3] ** 2)
        print(z[1])
        derivs=np.array([[z[2], -k*v**(N-1)*z[2] ,z[3], -g-k*v**(N-1)*z[3]]])
        return derivs
psolnE = RK4(f,z0,t,1,.1,x_0,vx_0,y_0,vy_0)```

Solution

  • The following statement always creates a numpy array of size 1, where derivs[0] being [z[2], -k*v**(N-1)*z[2] ,z[3], -g-k*v**(N-1)*z[3]].

    derivs=np.array([[z[2], -k*v**(N-1)*z[2] ,z[3], -g-k*v**(N-1)*z[3]]])
    

    Hence, as soon as your for loop in RK4 reaches 1, it throws an error. To resolve this, replace the above statement with-

    derivs=np.array([z[2], -k*v**(N-1)*z[2] ,z[3], -g-k*v**(N-1)*z[3]])