使用在每个时间步调用的函数求解微分方程

Question

方程组

def SEIRD_gov(y, t, beta_0, c_1, c_2, sigma, gamma, dr, ro, tg, g_1, g_2, ind):
    S, E, I, R, D = y
    
    dSdt = -beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2, ind) * S * I/N
    dEdt = beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2, ind) * I * S/N - sigma * E
    dIdt = sigma * E - (1 - dr) * gamma * I - dr * ro * I
    dRdt = (1 - dr) * gamma * I
    dDdt = dr * ro * I

    return dSdt, dEdt, dIdt, dRdt, dDdt

beta_gov - 这也是一个函数

def beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2, ind):
    beta_rez = beta_0 * gov(t, tg, g_1, g_2) * c_sig(t, c_1, c_2, ind)
    return beta_rez

但它也调用了两个函数

gov - 功能没问题

def gov(t, tg, g_1, g_2):
    if t > tg:
        alpha = 1 - g_1
    else:
        alpha = 1 - g_2
    return alpha

那个我不明白。

def c_sig(t, c_1, c_2, ind):
    sig = 1 / (1 + math.exp(c_1*(ind - c_2)))
    return sig

ind - DataSeries，一组数值，其中有很大一部分。当我调用主函数“SEIRD_gov”时，这些来自“ind”的值必须一次加一个来求解方程，然后将结果转移到微分方程组。

ind = df_region['self_isolation'].apply(lambda x: int(x))
ind = ind.values

可能在这里你需要添加一些类似循环的东西，但我不明白当一个函数被另一个函数调用时如何做到这一点。

算法应该如下：

beta_gov——调用两个函数，一个元素依次从ind加到c_sig，求解微分方程的returns个值。

以前，我只能传递此列表中的一项，这导致了错误的解决方案。

之前是Julia的模型，这里是部分带方程的代码。我只需要按已经存在的参数显示图形，但首先我需要在 python

中重写模型

y0 =  S0, E0, I0, R0, D0
ret = odeint(SEIRD_gov, y0, t, args=(beta_0, c_1, c_2, sigma, gamma, dr, ro, tg, g_1, g_2, ind))
S, E, I, R, D = ret.T

function SEIRD_gov!(du,u, p, t)
    S,E,I,R,D = u
    beta_0, c_1, c_2, sigma, gamma, dr, ro, tg, g_1, g_2 = p
    
    du[1] = -beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2) * S * I/N
    du[2] = beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2) * I * S/N - sigma * E
    du[3] = sigma * E - (1 - dr) * gamma * I - dr * ro * I
    du[4]= (1 - dr) * gamma * I 
    du[5] = dr * ro * I
end

function si(t)
    ind = convert(Int, round(t + 1))
    return data.self_isolation[ind]
end

c_lin(t, c_1, c_2) = 1 + c_1*(1 - c_2*si(t))
c_sig(t, c_1, c_2) = 1/(1 + exp(c_1*(si(t) - c_2)))

function gov(t, tg, g_1, g_2) 
    if t > tg
        alpha = 1 - g_1
    else
        alpha = 1 - g_2
    end
    alpha
end

beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2) = beta_0 * gov(t, tg, g_1, g_2)* c_sig(t, c_1, c_2)

结果图：

Answer 1

我猜你错过的是 ind 必须是一个函数，因为它是一个随时间变化的系数。

我假设 df_region 是来自 pandas 的 DataFrame。

因此，将代码中的 ind 定义修改为：

self_isolation_values = df_region['self_isolation'].to_numpy()
ind = lambda t: self_isolation_values[int(t)]

这使得ind成为t的函数，这样ind(t)就会将t转换为int和return相应的值您从 df_region['self_isolation'] 中的输入数据中获得的系数。这正是 julia 代码中函数 si(t) 的行为。

然后，在 c_sig 函数中，调用 ind(t) 函数

def c_sig(t, c_1, c_2, ind):
    sig = 1 / (1 + math.exp(c_1*(ind(t) - c_2)))
    return sig

使用在每个时间步调用的函数求解微分方程

Solving differential equations with a function called at each time step

python

function

julia

differential-equations