RK45 code w an error that says "index 50 is out of bounds for axis 0 with size 50"

Question

Here is my code:

import numpy as np
import time
import matplotlib as plt

def fun1(t,y):
    f=np.exp(-6*t)
    return (f)

def RK45Classic(h,f,t,yold):
    k1=h*f(t,yold)
    k2=h*f(t+h/2,yold+k1/2)
    k3=h*f(t+h/2,yold+k2/2)
    k4=h*f(t+h,yold+k3)
    ynew=yold+(1/6)*(k1+2*k2+2*k3+k4)
    return ynew

T_f=10
h=0.001
Nt=int(T_f/h)+1

t=np.linspace(0,T_f<Nt)
y=np.zeros(t.size)
f=np.zeros(t.size)

t0=time.time()
y[0]=1.
f[0]=fun1(t[0],y[0])

y[1]=RK45Classic(h,fun1,t[0],y[0])
f[1]=fun1(t[1],y[1])

y[2]=RK45Classic(h,fun1,t[1],y[1])
f[2]=fun1(t[2],y[2])

y[3]=RK45Classic(h,fun1,t[2],y[2])
f[3]=fun1(t[3],y[3])

for i in range(3,Nt):
    y[i+1]=y[i]+(h/24)*(55*f[i]-59*f[i-1]+37*f[i-2]-9*f[i-3])
    f[i+1]=fun1(t[i+1],y[i+1])

tfinal=time.time()-t0

plt.plot(t,y,t,np.exp(-6*t))

print(tfinal)

here is the error I get

IndexError                
Traceback (most recent call last)
<ipython-input-16-e315d11e30e6> in <module>
     35 
     36 for i in range(3,Nt):
---> 37     y[i+1]=y[i]+(h/24)*(55*f[i]-59*f[i-1]+37*f[i-2]-9*f[i-3])
     38     f[i+1]=fun1(t[i+1],y[i+1])
     39 

IndexError: index 50 is out of bounds for axis 0 with size 50

Some index in the error line is getting too large, probably the `[i+1]`. The upper bound of your `range(3,Nt)` is too large. If `y` has shape (50,), the index can't be 50 or larger. — hpaulj, Mar 07 '20 at 16:52

Lutz Lehmann · Answer 1 · 2020-03-07T17:36:29.480

What you implement in RK45Classic is the classical Heun-Kutta 4th order Runge-Kutta method, it is neither of Cash-Karp, Fehlberg nor Dormand-Prince embedded 45 methods.

A first measure would be to try to see why that error happens, so before starting the loop put

print(t.size, t[0], t[-1], y.shape)

This should show you that something in the construction of t went wrong.

The apparent error is

t=np.linspace(0,T_f<Nt)

correcting that typo should clean up that error. Effectively it creates a subdivision of [0,1] with a default size, apperently 50. Alternatively, why not use

t = np.arange(0, T_f+0.1*h, h)
Nt = len(t)-1

as that fits your input data more closely and serves the same purpose.

You might still get a similar error as there is no element with index Nt in the arrays t and y, the loop should iterate over range(3,len(t)-1), and if you want to put Nt there, its expression should give the same value, so that y[i+1] always addresses a valid array element.

Apart from that, the 4th order Adams-Bashford iteration itself looks fine.

RK45 code w an error that says "index 50 is out of bounds for axis 0 with size 50"

1 Answers1