python - 使用 ODEINT 或其他方法求解许多耦合微分方程组

Question

我正在尝试解决由 N=400 个神经元 组成的网络的动力学问题。

这意味着我有 400 个遵守以下规则的耦合方程:

i = 0,1,2...399
J(i,j) = some function of i and j (where j is a dummy variable)
I(i) = some function of i
dr(i,t)/dt = -r(i,t) + sum over j from 0 to 399[J(i,j)*r(j)] + I(i)

我该如何解决？

我知道对于 3 首颂歌的系统。我定义了 3 个 odes 和初始条件，然后应用 odeint。在这种情况下是否有更好的执行方式？

到目前为止，我尝试了以下代码(它并不好，因为它进入了无限循环):

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()

这是脚本生成的情节: