Python Odeint дает странные результаты - PullRequest
Новая
       70

Python Odeint дает странные результаты

0 голосов
/ 11 сентября 2018

Я пытаюсь понять, как работает scipy.odeint, но у меня есть некоторые проблемы. Для экзамена: 

import numpy as np
from scipy.integrate import  odeint
from scipy.integrate import  odeint

def Diffeq(v,t, lam, gam,c, a):
    vdot = v
    for i in range(0,len(v)):
        if i == 0:
            vdot[0] = c[1]*v[1]- lam*a[0]*v[0]
        elif i == (len(v)-1):
            vdot[i] =  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]
        else:
            vdot[i]=  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]+ c[i+1]*v[i+1]
    print vdot
    return vdot



incond=np.array([0]*900)
incond[1] =1
t = np.linspace(0.0, 2, 1000)
ak = [2*i for i in range(0,900)]
lamma =2
gamma =1
c=[i*gamma for i in range(0,900)]

y = odeint(Diffeq, incond, t, args=(lamma,gamma,c,ak) )

Этот код должен вычислять систему дифференциальных уравнений вида: enter image description here

, где

xdot_0 = - (a_0 + c_0) * x_0 (t) + c_1 * x_1 (t)
xdot_899 = a_898 * x_898 (т) - (a_899 + c_899) * x_900 (т)

с начальным условием x (0) = (0,1,0 ... 0) когда я пытаюсь проанализировать результаты, я замечаю, что мои функции взрываются до + бесконечности. Если я поиграю с константами ak, lamma и gamma, я могу получить результаты, похожие на

x_k (t) = [0, -21,21, 0, 0, ..., 0]

при каждой операции. Поэтому я думаю, что допустил ошибку в своем коде, но не вижу, где.

1 Ответ

0 голосов
/ 11 сентября 2018

В Python при выполнении строки 

    vdot = v

vdot не является копией v. Два имени теперь относятся к одному и тому же объекту. Поэтому, когда вы изменяете vdot в функции Diffeq, вы также изменяете входной аргумент. Если вы измените vdot[0], а затем попытаетесь использовать v[0], вы фактически получите измененное значение vdot[0], поэтому ваши вычисления неверны.

Измените эту строку, скажем, 

    vdot = np.empty_like(v)

Когда я это делаю (и я удаляю оператор print, чтобы функция завершалась в разумные сроки), odeint успешно возвращается. Вот график компонентов решения:

plot

...