问题

在我的代码中，我希望gamma的值从0.4逐渐增加到0.8，可以说是线性增加，而且增加的步幅很小。我希望随着时间t的增加，gamma的值也随之增加，所以当时间序列结束时，gamma的值是0.8。在这个常微分方程系统中是否可以实现这个目标呢？

# 总人口，N。
N = 1
# 初始感染和康复人数，I0和R0。
I0, R0 = 0.001, 0
# 其他人，S0，在最初时刻都容易受到感染。
U0 = N - I0 - R0
J0 = I0
Lf0, Ls0 = 0, 0
# 接触率，beta，和平均康复率，gamma，（以1/天为单位）。
beta, gamma = 8, 0.4
int_gamma = np.linspace(0.4, 0.8, 1000+1)
mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
t = np.linspace(0, 1000, 1000+1)

# SIR模型的微分方程。
def deriv(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
    dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI


# 在时间网格t上积分SIR方程。
solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
U, Lf, Ls, I, R, cInc = solve.T

import numpy as np
import matplotlib.pyplot as plt



# Total population, N.
N = 1
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 0.001, 0
# Everyone else, S0, is susceptible to infection initially.
U0 = N - I0 - R0
J0 = I0
Lf0, Ls0 = 0, 0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 8, 0.4
int_gamma = np.linspace(0.4, 0.8, 1000+1)
mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
t = np.linspace(0, 1000, 1000+1)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
    dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI


# Integrate the SIR equations over the time grid, t.
solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
U, Lf, Ls, I, R, cInc = solve.T

答案1

得分: 1

There is np.geomspace which may help.

import numpy as np

arr = np.geomspace(0.4, 0.8, num=10)
arr
# array([0.4       , 0.4320239 , 0.46661162, 0.50396842, 0.544316,
#        0.5878938 , 0.63496042, 0.68579519, 0.74069977, 0.8       ])
np.diff(arr)
# array([0.0320239 , 0.03458772, 0.0373568 , 0.04034758, 0.0435778 ,
#        0.04706662, 0.05083477, 0.05490458, 0.05930023])

Or compound two linspaces

x = 1.5
arr = np.linspace(0.4, 0.8/x, num=10) * np.linspace(1, x, num=10)
arr
# array([0.4       , 0.43786008, 0.47736626, 0.51851852, 0.56131687,
#        0.60576132, 0.65185185, 0.69958848, 0.74897119, 0.8       ])

np.diff(arr)
# array([0.03786008, 0.03950617, 0.04115226, 0.04279835, 0.04444444,
#        0.04609053, 0.04773663, 0.04938272, 0.05102881])

Different values of x change the acceleration in the result.

<details>
<summary>英文:</summary>

There is `np.geomspace` which may help.

    import numpy as np

    arr = np.geomspace( 0.4, 0.8, num = 10 )
    arr
    # array([0.4       , 0.4320239 , 0.46661162, 0.50396842, 0.544316, 
    #        0.5878938 , 0.63496042, 0.68579519, 0.74069977, 0.8       ])
    np.diff( arr )
    # array([0.0320239 , 0.03458772, 0.0373568 , 0.04034758, 0.0435778 ,
    #        0.04706662, 0.05083477, 0.05490458, 0.05930023])

Or compound two linspaces
    
    x = 1.5
    arr = np.linspace( 0.4, 0.8/x, num = 10 ) * np.linspace( 1, x, num = 10 )
    arr
    # array([0.4       , 0.43786008, 0.47736626, 0.51851852, 0.56131687,
    #        0.60576132, 0.65185185, 0.69958848, 0.74897119, 0.8       ])

    np.diff( arr )  
    # array([0.03786008, 0.03950617, 0.04115226, 0.04279835, 0.04444444,
    #        0.04609053, 0.04773663, 0.04938272, 0.05102881])
   
Different values of x change the acceleration in the result.

    

</details>