Solving a system of many coupled differential equations using ODEINT or something else

Question

I am trying to solve the dynamics of a network composed of N=400 neurons.

That means I have 400 coupled equations that obey the following rules:

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)

How do I solve?

I know that for a system of 3 odes. I defined the 3 odes and the initial conditions and then apply odeint. Is there a better way to perform in this case?

So far I tried the following code (it isn't good since it enters an infinite loop):

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)

In the future, it would be helpful if the code that you show in the question can be run without being edited (see "How to create a Minimal, Complete, and Verifiable example"). Currently the code is missing import statements and has some undefined variables. — Warren Weckesser
– Warren Weckesser, Commented May 5, 2019 at 18:58

Warren Weckesser · Accepted Answer · 2019-05-05 19:07:46Z

After making what I think are the obvious corrections and additions to your code, I was able to run it. It was not actually in an infinite loop, it was just very slow. You can greatly improve the performance by "vectorizing" your calculations as much as possible. This allows the loops to be computed in C code rather than Python. A hint that there is room for a lot of improvement is in the expression sum over j from 0 to 399[J(i,j)*r(j)]. That is another way of expressing the product of a matrix J and a vector r. So we should really have something like J @ r in the code, and not all those explicit Python loops.

After some more tweaking, here's a modified version of your code. It is significantly faster than the original. I also reorganized a bit, and added a plot.

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

Here's the plot generated by the script:

hello warren , i had a follow up question that you might be able to answer . i posted it here stackoverflow.com/questions/56039218/…

Collectives™ on Stack Overflow

Solving a system of many coupled differential equations using ODEINT or something else

1 Answer 1

1 Comment

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Related