Решение системы многих связанных дифференциальных уравнений с использованием ODEINT или чего-то еще

Question

Я пытаюсь решить динамику сети, состоящей из N = 400 нейронов .

Это означает, что у меня есть 400 связанных уравнений, которые подчиняются следующим правилам:

i = 0,1,2 ... 399
J (i, j) = некоторая функция от i и j (где j - фиктивная переменная)
Я (я) = некоторая функция от меня
dr (i, t) / dt = -r (i, t) + сумма по j от 0 до 399 [J (i, j) * r (j)] + I (i)

Как мне решить?

Я знаю, что для системы из 3 од. Я определил 3 оды и начальные условия, а затем применил одеинт. Есть ли лучший способ выполнить в этом случае?

До сих пор я пробовал следующий код (он не очень хорош, поскольку он входит в бесконечный цикл):

N=400
t=np.linspace(0,20,1000)
J0=0.5
J1=2.5
I0=0.5
I1=0.001
i=np.arange(0,400,1)
theta=(2*np.pi)*i/N
I=I0+I1*cos(theta)
r=np.zeros(400)
x0 = [np.random.rand() for ii in i]

def threshold(y):
    if y>0:
        return y
    else:
        return 0

def vectors(x,t):
    for ii in i:
        r[ii]=x[ii]
    for ii in i:
        drdt[ii] = -r[ii] + threshold(I[ii]+sum(r[iii]*(J0+J1*cos(theta[ii]-theta[iii]))/N for iii in i))
    return drdt

x=odeint(vectors,x0,t)

Answer 1

После внесения очевидных исправлений и дополнений в ваш код я смог его запустить.Это не было на самом деле в бесконечном цикле, это было просто очень медленно.Вы можете значительно улучшить производительность, максимально «векторизовав» свои расчеты.Это позволяет вычислять циклы в C-коде, а не в Python.Намек на то, что есть место для значительных улучшений, содержится в выражении sum over j from 0 to 399[J(i,j)*r(j)].Это еще один способ выразить произведение матрицы J и вектора r.Таким образом, в коде должно быть что-то вроде J @ r, а не все эти явные циклы Python.

После некоторой дополнительной настройки, вот модифицированная версия вашего кода.Это значительно быстрее, чем оригинал.Я также немного реорганизовал и добавил график.

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


def computeIJ(N):
    i = np.arange(N)
    theta = (2*np.pi)*i/N

    I0 = 0.5
    I1 = 0.001
    I = I0 + I1*np.cos(theta)

    J0 = 0.5
    J1 = 2.5
    delta_theta = np.subtract.outer(theta, theta)
    J = J0 + J1*np.cos(delta_theta)
    return I, J / N


def vectors2(r, t, I, J):
    s = J @ r
    drdt = -r + np.maximum(I + s, 0)
    return drdt


N = 400

I, J = computeIJ(N)

np.random.seed(123)
r0 = np.random.rand(N)
t = np.linspace(0, 20, 1000)

r = odeint(vectors2, r0, t, args=(I, J))

for i in [0, 100, 200, 300, 399]:
    plt.plot(t, r[:, i], label='i = %d' % i)

plt.xlabel('t')
plt.legend(shadow=True)
plt.grid()

plt.show()

Вот график, сгенерированный скриптом: