I would like to solve a system of network-based differential equations using python 3.6

Question

I would like to solve a system of network-based differential equations using python 3.6. The system of equations is as follows:

dx_i/dt = omega_i  - epsilon_i * y_i * x_i - mu_i * x_i,  

dy_i/dt = epsilon_i * y_i * x_i - zeta_i * y_i - rho_i * y_i * z_i,

dv_i/dt = c_i * y_i - gamma_i * v_i + \sum_{i neq j} beta_{ij} * v_i * x_i,

dz_i/dt = k_i * y_i * z_i - delta_i * z_i,

where beta_{ij} = beta (1 - sigma_{ij}) * exp(- alpha|i-j|) 

i = 1,2,3,...,N

I have written the code below in an attempt to solve the system of differential equations in a spatial network.

from jitcode import jitcode, y
import numpy as np
import sympy

#import symengine
import matplotlib.pyplot as plt
#from scipy.integrate import ode
from numpy.random import uniform

n = 10
alpha = 0.05
#beta = 0.1

beta = uniform(0.01,3.0,n)
beta.sort()
mu = uniform(0.01,3.0,n)
mu.sort()
epsilon = uniform(0.01,3.0,n)
epsilon.sort()
pi = uniform(0.01,3.0,n)
pi.sort()
gamma = uniform(0.01,3.0,n)
gamma.sort()
omega = uniform(0.01,3.0,n)
omega.sort()
zeta = uniform(0.01,3.0,n)
zeta.sort()
rho = uniform(0.01,3.0,n)
rho.sort()
k = uniform(0.01,3.0,n)
k.sort()
c = uniform(0.01,3.0,n)
c.sort()


# Knonecker delta

M = np.einsum('ij-> ij',np.eye(n,n))

print(M)

# Adjacency matrix

A = beta * M * sympy.exp(-alpha) 

print(A)


def f():
    for i in range(n):

    coupling_sum = sum(y(i+2) * y(i) for j in range(n) if A[i, j]
    )
    yield omega[i] - epsilon[i] * y(i+2) * y(i) - mu[i] * y(i)
    yield epsilon[i] * y(i+2) * y(i) - zeta[i] * y(i+1) - rho[i] * y(i+1)* y(i+3)
    yield c[i] * y(i+1) - gamma[i] * y(i+2) + coupling_sum
    yield k[i]* y(i+1) * y(i+3) - delta[i] *y(i+3)

#integrate
#---------------

initial_state = np.random.random(n)

ODE = jitcode(f,n=n)
ODE.set_integrator("dopri5", atol=1e-6,rtol=0)
initial = np.linspace(0.1,0.4,num=n)

ODE.set_initial_value(initial_state,time=0.0)
#data structure: x[0], w[0], v[0], z[0], ..., x[n], w[n], v[n], z[n]
data = []
data = np.vstack(ODE.integrate(T) for T in range(0, 100, 0.1))
print(data)
fig = plt.figure()

I am expecting to get solutions for the four differential equations and some simulations to represent the equations. The error message that Iam getting says "RuntimeError: Not Implemented"

Answer 1

Your translation of the model into code has index problems. You addressed this once in the sum calculation, but then never again. To help with that, define helper functions

def X(i): return y(4*i)
def Y(i): return y(4*i+1)
def V(i): return y(4*i+2)
def Z(i): return y(4*i+3)

Then you can express the generator of the symbolic right sides as

def f():
    for i in range(n):

        coupling_sum = V(i) * sum(beta[i,j]*X(j) for j in range(n) if j!=i )
        yield omega[i] - epsilon[i] * Y(i) * X(i) - mu[i] * X(i)
        yield epsilon[i] * Y(i) * X(i) - zeta[i] * Y(i) - rho[i] * Y(i)*Z(i)
        yield c[i] * Y(i) - gamma[i] * V(i) + coupling_sum
        yield k[i] * Y(i) * Z(i) - delta[i] * Z(i)

with an appropriate definition of beta[i,j]. The code you wrote is strange and does not match the formula. In the formula, beta is first a constant and then a matrix. In the code, beta is an array. This is quite incompatible.

In the call to compile the function, you also should give the correct dimension of the state space, you have n components of x,y,v,z, making 4*n components in total.

initial_state = np.random.random(4*n)
ODE = jitcode(f,n=4*n)

With the correct dimension of the state space arguments, the integration calls probably will go through.