-2

我想使用 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”

4

1 回答 1

1

您将模型转换为代码时存在索引问题。您在总和计算中解决了这个问题,但以后再也没有解决过。为了解决这个问题,定义辅助函数

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)

然后您可以将符号右侧的生成器表示为

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)

有适当的定义beta[i,j]。你写的代码很奇怪,和公式不符。在公式中,beta首先是一个常数,然后是一个矩阵。在代码中,beta是一个数组。这是非常不兼容的。

在编译函数的调用中,您还应该给出状态空间的正确维度,您有n的组件,总共x,y,v,z制作4*n组件。

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

使用状态空间参数的正确维度,集成调用可能会通过。

于 2019-08-13T06:19:54.983 回答