Backward bifurcation

3121 days ago by medlock

var('Pi_0 N alpha beta mu nu epsilon S I V N s i v') numbers = [Pi_0 * N - beta * S * I / N - alpha * S - mu * S, alpha * S - (1 - epsilon) * beta * V * I / N - mu * V, beta * (S + (1 - epsilon) * V) * I / N - (mu + nu) * I, (Pi_0 - mu) * N - nu * I] show(numbers) 
       
X = (S, V, I) proportions = [expand((numbers[k] / N - numbers[-1] / N * X[k] / N).subs(S = s * N, V = v * N, I = i * N)) for k in xrange(len(X))] show(proportions) 
       
x = (s, v, i) J = jacobian(proportions, x) show(J) 
       
DFE = {s: Pi_0 / (alpha + Pi_0), v: alpha / (alpha + Pi_0), i: 0} J0 = J.subs(DFE) show(J0) 
       
sigma0 = J0.eigenvalues() show(sigma0) 
       
beta_C = (solve(sigma0[0], beta)[0]).rhs() show(beta_C) 
       
ev = J0.eigenmatrix_left()[1][0, :] ew = J0.eigenmatrix_right()[1][:, 0] n = (ev * ew)[0] ew = vector([x / n for x in ew]).column() 
       
x = (s, v, i) b = sum([ sum([ ev[0, k0] * ew[k1, 0] * diff(proportions[k0], x[k1], beta).subs(DFE).subs(beta = beta_C) for k1 in xrange(len(x))]) for k0 in xrange(len(x))]) show(b) 
       
x = (s, v, i) a = 1/2 * sum([ sum([ sum([ (ev[0, k0] * ew[k1, 0] * ew[k2, 0] * diff(proportions[k0], x[k1], x[k2]).subs(DFE).subs(beta = beta_C)) for k2 in range(len(vars))]) for k1 in range(len(vars))]) for k0 in range(len(vars))]) show(a)