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)
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)