# Solving systems of nonlinear ODEs in epidemiology

I have a 9 systems of nonlinear ODEs to solve. I want to determine the endemic equilibrium points.

How do I go about it?

I tried manual calculation but it becomes cumbersome. Can software be used for this purpose?

• Do you need an analytical solution or will it suffice to have a numerical solution for a particular choice of parameters in your ODE? – Paul Aug 7 '17 at 18:29
• Yes, I will prefer an analytical solution.. N which can you help with? – Kareem Adams Aug 7 '17 at 18:34
• With analytical solutions, I can then see the behaviour of the system as parameters are being altered.. – Kareem Adams Aug 7 '17 at 18:36
• Use Mathematica or something like SymEngine in Julia. There's a good change there is no analytical solution though, so numerical solutions may be the only option. – Chris Rackauckas Aug 7 '17 at 20:24
• Endemic equilibrium(i.stack.imgur.com/DQq8v.jpg Here s d equation)](i.stack.imgur.com/DQq8v.jpg Here s d equation) – Kareem Adams Aug 8 '17 at 11:12

If you want to find the equilibrium points you need to write your system as a first-order one

$$\frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x})$$

and solve the non-linear system of equations

$$\mathbf{F}(\mathbf{x})=\mathbf{0}\ ,$$

that might not be solvable analytically. After that, you need to find the eigenvalues of the Jacobian matrix around each equilibrium point $$\mathbf{x}^*$$, $$\mathbf{J}(\mathbf{x}^*)$$. You can use a Computer Algebra System (CAS) for the job. SymPy, Maxima, Mathematica and Maple are some examples.

An example of this, using SymPy, for the Lotka-Volterra model is the following

from sympy import symbols, Matrix, solve, diff
x, y = symbols("x y")
alpha, beta, gamma, delta = symbols("alpha beta gamma delta")
dXdt = Matrix([alpha*x - beta*x*y, delta*x*y -gamma*y])
eqpts = solve(dXdt, [x, y])
J = Matrix(2, 2, lambda i, j: diff(dXdt[i], [x, y][j]))
vals1 = J.subs({x: eqpts[0][0], y: eqpts[0][1]}).eigenvals()
vals2 = J.subs({x: eqpts[1][0], y: eqpts[1][1]}).eigenvals()
print(vals1)
print(vals2)


with output

{-gamma: 1, alpha: 1}
{sqrt(-alpha*gamma): 1, -sqrt(-alpha*gamma): 1}


Here you should take into account two things:

1. The system of non-linear equations might not have an analytical solution, and you might need to solve it numerically (with Newton method, for example).
2. Even if you have the analytical solution(s) for your system, the eigensystem might not be solvable analytically.