Код RK45 с ошибкой, которая говорит: «Индекс 50 выходит за пределы для оси 0 с размером 50» - PullRequest
Новая
       64

Код RK45 с ошибкой, которая говорит: «Индекс 50 выходит за пределы для оси 0 с размером 50»

0 голосов
/ 07 марта 2020

Вот мой код: 

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

1 Ответ

0 голосов
/ 07 марта 2020

То, что вы реализуете в RK45Classic, - это классический метод Рунге-Кутты 4-го порядка Хеун-Кутты, в котором нет ни методов, встроенных Ка sh -Карпом, Фельбергом или Дорманд-Принсом.

A Первой мерой было бы попытаться выяснить, почему происходит эта ошибка, поэтому перед запуском l oop put 

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, l oop должен перебрать range(3,len(t)-1), и если вы хотите поставить Nt там его выражение должно давать одно и то же значение, так что y[i+1] всегда обращается к допустимому элементу массива.

Кроме того, сама итерация Адамса-Башфорда 4-го порядка выглядит нормально.

...