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

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

但它也调用了两个函数

1. gov - 功能没问题

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

1. 我不明白那个。

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添加一个元素，并返回求解微分方程的值。

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

之前它是 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 并从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