Решаем комплексную систему оду с одеинтом Python - PullRequest
Новая
       93

Решаем комплексную систему оду с одеинтом Python

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

Я пытаюсь решить эту систему:

P ̇ (t) = [γ - δ - α - u (t)] P (t) + βQ (t) (1) Q ̇ ( t) = αP (t) - (λ + β) Q (t) (2)

, где u - функция с уровнями (u (t) = lesu (i), где i таково, что t ∈ [(i - 1) T / N, i T / N [от lesu, который является вектором N последовательных управляющих действий):

https://i.stack.imgur.com/c5lF7.png

Я определил функцию для вычисления u, системы и попытался решить эту систему с помощью OdeInt, как показано ниже 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

alph=5.63
bet=0.48
gamm=1.46
delt=0.16
T=30
m=0.001
lambd=0.16
rh=0.35
lesu=[0.5 ,0.8 ,0.3 ,1.9 ,0.54]
compt=0
time=np.linspace(0,30)

def udet(t,lesu):
    j=0 
    u=0
    N=len(lesu)
    lest=[]
    for i in range (N):
        lest.append(i*(T/N))
    j=lest.index(max(x for x in lest if x<=t))
    u=lesu[j]
    return u

print(udet(T/3,lesu))

def my_ode(der,t):
    global compt
    u=udet(time[compt],lesu)
    print(time[compt])
    compt=compt+1
    if compt<len(time):
        print(compt)
        res1=(gamm-delt-alph-u)*der[0]+bet*der[1]
        res2=alph*der[0]-(lambd+bet)*der[1]
        res=[res1,res2]
        return res


P0=0.5
Q0=0.5
res0=[P0,Q0]

y = odeint(my_ode,res0,time)

plt.plot(time,y[:,0])
plt.plot(time,y[:,1])

Однако это не работает из-за вычисления u в системе , У меня есть эта ошибка: 

    u=udet(time[compt],lesu)

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

Я действительно не понимаю, почему у меня есть эта ошибка? Должен ли я использовать другой метод решения?

Заранее спасибо

...