python - 使用具有不同参数的 solve_ivp

Question

我正在尝试求解以下简单的线性方程组:

x'(t) = A_eps(t) x(t)

其中 x 是一个 n 维向量，A_eps(t) 是一个随时间变化的矩阵。

以下是我在帖子之后所做的尝试:

scipy ode update set_f_params inside function set as set_solout

首先我将右侧定义为一个函数:

from functools import partial
from scipy.integrate import solve_ivp

def Aeps_x(t, x, enviroment, Q, F, H):

    '''
    inputs: 
        - t = time.
        - x[j] = # of cell with phenotype j in an enviroment. 
        - enviroment = number of enviroment. There are E in total.  
        - Q = ExE transition matrix of a Markov Chain.  
        - F = ExP matrix with growth rate of phenotypes for each enviroment. 
        - H = PxP matrix with the switching rates from phenotype i to j.
    '''
    E = Q.shape[0]
    P = F.shape[1]

    A = np.diag(F[enviroment]) - H

    dx = A.dot(x)

    return(dx)

然后，我只是为 r.h.s. 设置了一些参数:

EMat = np.array([[0, 1, 1, 1], 
             [1, 0, 1, 1],
             [1, 1, 0, 1],
             [1, 1, 1, 0]])
E0 = EMat.shape[0]

row_sums = EMat.sum(axis=1).reshape(E0, 1)
Q = EMat / row_sums

F = linalg.toeplitz([1, 0.1, 0.1, 0.1]) # only one strong phenotype for each 
enviroment
F = F[0:E0, ] # only 4 enviroments

P0 = F.shape[1]
H = np.diag([0.5]*P0)

为了设置求解器，我做了:

sol_param = solve_ivp(
    partial(Aeps_x, enviroment = 2, Q=Q, F=F, H=H), (0, 4, 20), x0, dense_output=True)

我想写这样的东西:

sol_param = solve_ivp(
    partial(Aeps_x, enviroment = next_enviroment(current_env, Q), 
    Q=Q, F=F, H=H), (0, 4, 20), x0, dense_output=True)

其中 next_enviroment(current_env, Q) 是:

def next_enviroment(current_env, Q):
    '''
    Given an current state, computes the next state in a markov chain with transition matrix Q. 
    '''
    sample = np.random.multinomial(1, Q[intitial_env], size = 1)
    next_env = np.argmax(sample)
    return(next_env)

是一个函数，它获取当前状态并根据给定规则选择一个随机的新状态。我的问题有两个:

• 首先，我不知道如何在求解器中读取当前环境。
• 其次，给定当前环境，如何将其传递到函数中。

感谢您的帮助。

Answer 1

我找到了答案。这是代码:

EMat = np.array([[10, 0.1], 
                 [0.1, 10]])

# number of enviroments and phenotypes
E = 2
P = 2

row_sums = EMat.sum(axis=1).reshape(E, 1)
Q = EMat / row_sums

H = np.array([[0, 0.05],
              [1e-6, 0]])

F = np.array([[2, -0.05], 
              [-7, -0.05]]) 


import scipy

N = 1000
tfinal = 25
t = np.linspace(start=0, stop=tfinal, num=N)
t0 = 0
x0 = [20, 20]
e0 = 0


solver = scipy.integrate.ode(Aeps_x).set_integrator('dopri5', nsteps=100)
solver.set_initial_value(x0, t0).set_f_params(e0, Q, F, H)

sol = np.zeros((N, E))
sol[0] = x0
Enviroments = np.zeros(N)

k = 1
while solver.successful() and solver.t < tfinal:
    solver.integrate(t[k])
    sol[k] = solver.y
    k += 1

    Enviroments[k-1] = e0

    e0 =  next_enviroment(e0, Q=Q)
    solver.set_f_params(e0, Q, F, H)

这是模拟的 pplot:

Evolution of two phenotypes in an stochastic random enviroment with two states. State j favors phenotype j