In my code below, I wish to have the value of gamma increase from 0.4, to 0.8, lets say linearly at first, in small increments. I want this to go up as time t goes up, so by the time the time series is over, gamma is 0.8. Is it possible to do this within this ODE system?
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
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.