0

我想为在 Python 3.6 中连接到网络中的个人绘制微分方程组的图表。方程组如下：

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

dy_i/dt = r_i * y_i - gamma * x_i * y_i + \sum_{j!=1} A_{ji} *y_i


x_i(t) is the antibody response in the i-th individual
y_i(t) is the viral charge in that individual where i = 1,....,n
omega_i is the rate of production and/or transport of antibodies
mu_i is the death rate of antibodies
epsilon_i is the rate of production of antibodies induced by a unit viral 
population
r_i is the intrinsic growth rate of viral population
gamma_i is the rate of destruction of viruses by a unit antibody population
A_{ji} is the ji-th of a matrix representing the strength of transmission from j to i

我已经为连接在网络中的n 个个体编写了对病毒入侵的免疫反应的代码。

该模型代表了一个耦合方程系统，表示连接个体网络中的抗体和病毒种群。

from jitcode import jitcode, y
import numpy as np
import sympy
import matplotlib.pyplot as plt
from numpy.random import uniform
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D


n = 5
alpha = 0.05


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()
gamma = uniform(0.01,3.0,n)
gamma.sort()
omega = uniform(0.01,3.0,n)
omega.sort()
r = uniform(0.01,3.0,n)
r.sort()


# Knonecker delta

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

print(M)

# The transmission matrix A whose elements represent the strength of 
# transmission from j to i depending of spatial factors. 

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



def X(i): return y(2*i)
def Y(i): return y(2*i+1)


def f():
    for i in range(n):
    
        coupling_sum = sum(A[i,j]*Y(j) for j in range(n) if j!=i ) 
   
        yield omega[i] - mu[i] * X(i) + epsilon[i] * X(i) * Y(i)
        yield r[i] * Y(i) - gamma[i] * X(i) * Y(i) + coupling_sum
    
    

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

t = np.linspace(0, 100)
T = np.row_stack([t, t])

initial_state = np.random.random(2*n)

ODE = jitcode(f, n=2*n)
ODE.set_integrator("dopri5", atol=1e-6,rtol=0)
ODE.generate_f_C(simplify=False, do_cse=False, chunk_size=150)
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(10, 100, 10))
print(data)


# Plotting the graphs
#-----------------------

plt.show()
plt.savefig('tmp.pdf'); plt.savefig('tmp.png')
plt.title("The Immunoepidemiological model")
plt.plot(t, f)
plt.xlabel('t')
plt.ylabel('f')
fig = plt.figure()

我期望得到随着时间 t 的抗体和病毒种群的图表。但是，我收到以下错误消息。

ValueError：x 和 y 必须具有相同的第一维，但具有 (50,) 和 (1,) 形状

4

1 回答 1

1

您的大部分代码都很好，除非您没有绘制正确的东西。f是对您的微分方程的抽象描述，而不是您的解决方案。

您的集成数据包含在data变量中，其形状为 (timepoints, system_variables)。现在，例如，如果您想将抗体反应和病毒种群绘制为第一个个体的时间函数，您的绘图命令应该是：plt.plot(t, data[:,0:2]).

此外，您没有将t其用作集成的实际时间。按照 Jitcode文档，它应该是这样的：

t = np.arrange(0,100,0.1)
data = np.vstack(ODE.integrate(time) for time in t)
于 2019-08-16T18:32:58.467 回答