Backward bifurcation -- Sage

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)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{I S \beta}{N} + N \Pi_{0} - S \alpha - S \mu, \frac{{\left(\epsilon - 1\right)} I V \beta}{N} + S \alpha - V \mu, -{\left(\mu + \nu\right)} I - \frac{{\left({\left(\epsilon - 1\right)} V - S\right)} I \beta}{N}, {\left(\Pi_{0} - \mu\right)} N - I \nu\right]$

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)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\beta i s + i \nu s - \Pi_{0} s - \alpha s + \Pi_{0}, \beta \epsilon i v - \beta i v + i \nu v - \Pi_{0} v + \alpha s, -\beta \epsilon i v + \beta i s + \beta i v + i^{2} \nu - \Pi_{0} i - i \nu\right]$

x = (s, v, i) J = jacobian(proportions, x) show(J)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -\beta i + i \nu - \Pi_{0} - \alpha & 0 & -\beta s + \nu s \\ \alpha & \beta \epsilon i - \beta i + i \nu - \Pi_{0} & \beta \epsilon v - \beta v + \nu v \\ \beta i & -\beta \epsilon i + \beta i & -\beta \epsilon v + \beta s + \beta v + 2 \, i \nu - \Pi_{0} - \nu \end{array}\right)$

DFE = {s: Pi_0 / (alpha + Pi_0), v: alpha / (alpha + Pi_0), i: 0} J0 = J.subs(DFE) show(J0)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -\Pi_{0} - \alpha & 0 & -\frac{\Pi_{0} \beta}{\Pi_{0} + \alpha} + \frac{\Pi_{0} \nu}{\Pi_{0} + \alpha} \\ \alpha & -\Pi_{0} & \frac{\alpha \beta \epsilon}{\Pi_{0} + \alpha} - \frac{\alpha \beta}{\Pi_{0} + \alpha} + \frac{\alpha \nu}{\Pi_{0} + \alpha} \\ 0 & 0 & -\frac{\alpha \beta \epsilon}{\Pi_{0} + \alpha} + \frac{\Pi_{0} \beta}{\Pi_{0} + \alpha} + \frac{\alpha \beta}{\Pi_{0} + \alpha} - \Pi_{0} - \nu \end{array}\right)$

sigma0 = J0.eigenvalues() show(sigma0)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{\alpha \beta \epsilon + {\left(\alpha - \beta + \nu\right)} \Pi_{0} + \Pi_{0}^{2} - \alpha \beta + \alpha \nu}{\Pi_{0} + \alpha}, -\Pi_{0} - \alpha, -\Pi_{0}\right]$

beta_C = (solve(sigma0[0], beta)[0]).rhs() show(beta_C)

 $\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\alpha + \nu\right)} \Pi_{0} + \Pi_{0}^{2} + \alpha \nu}{\alpha \epsilon - \Pi_{0} - \alpha}$

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)

 $\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\alpha \epsilon}{\Pi_{0} + \alpha} + \frac{\Pi_{0}}{\Pi_{0} + \alpha} + \frac{\alpha}{\Pi_{0} + \alpha}$

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)

 $\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(\beta - \nu\right)} {\left({\left(\alpha + \nu\right)} \Pi_{0} + \Pi_{0}^{2} + \alpha \nu\right)} \Pi_{0}}{{\left(\alpha \epsilon - \Pi_{0} - \alpha\right)} {\left(\alpha \beta \epsilon - {\left(\alpha + \beta - \nu\right)} \Pi_{0} - \alpha^{2} - \alpha \beta + \alpha \nu\right)}} - \frac{{\left(\frac{{\left({\left(\alpha + \nu\right)} \Pi_{0} + \Pi_{0}^{2} + \alpha \nu\right)} \epsilon}{\alpha \epsilon - \Pi_{0} - \alpha} - \frac{{\left(\alpha + \nu\right)} \Pi_{0} + \Pi_{0}^{2} + \alpha \nu}{\alpha \epsilon - \Pi_{0} - \alpha}\right)} {\left(\alpha \beta \epsilon - \Pi_{0} \alpha - \alpha^{2} - \alpha \beta + \alpha \nu\right)}}{\alpha \beta \epsilon - {\left(\alpha + \beta - \nu\right)} \Pi_{0} - \alpha^{2} - \alpha \beta + \alpha \nu} + \nu$

Backward bifurcation

3121 days ago by medlock