RK45代码w表示“索引50超出大小为50的轴0的界限”的错误

Question

这是我的代码：

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)

这是我得到的错误

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

Answer 1

在{}中实现的是经典的Heun-Kutta四阶Runge-Kutta方法，它既不是Cash-Karp方法，也不是Fehlberg方法，也不是休眠Prince方法

第一项措施是尝试了解错误发生的原因，因此在开始循环前，请将

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

这应该表明t的构造中出现了一些错误

明显的错误是

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

纠正那个打字错误应该可以消除那个错误。实际上，它创建了一个[0,1]的细分，默认大小为50。或者，为什么不使用

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

因为这更适合您的输入数据，并且服务于相同的目的

您可能仍然会得到类似的错误，因为在数组t和y中没有索引为Nt的元素，循环应该在range(3,len(t)-1)上迭代，如果您想把Nt放在那里，它的表达式应该给出相同的值，以便y[i+1]始终寻址一个有效的数组元素

除此之外，四阶Adams-Bashford迭代本身看起来不错