我想使用 python 3.6 求解一个基于网络的微分方程系统。方程组如下:
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
我编写了下面的代码,试图解决空间网络中的微分方程组。
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()
我期待得到四个微分方程的解和一些模拟来表示方程。我收到的错误消息是“RuntimeError: Not Implemented”