Software to find the equilibria of a set of differential equations?

Question

Context for this question: I'm in mathematical biology, and while I've got a decent grasp of the concepts, and am a decent hand with simulation and numerical approaches, my analytical abilities and general skill at manipulating equations has always been rough. Probably the root of deciding I was "bad at math" back in the day, though that's neither here nor there.

What I'm trying to do is find the equilibriums for a set of differential equations, similar to the following:

$$\begin{align} \frac{{dS}}{{dt\strut}} =& \mu - (\beta I + \mu )S\\ \frac{{dE}}{{dt\strut}} =& \beta SI - (\mu + \sigma )E\\ \frac{{dI}}{{dt\strut}} =& \sigma E - (\mu + \gamma )I\\ \frac{{dR}}{{dt\strut}} =& \gamma I - \mu R \end{align}$$

I can do this one, though admittedly it takes me quite some time, but more complex systems involve a fair number of hours going down dead ends that really shouldn't be, basic errors etc.

Any chance any of the myriad math software packages out there have a way to simplify the computation/manipulation steps of finding the equilibria for a system like this? Any general tips if there aren't? You've got Matlab, Mathematica, and really any open source package you please to work with.

Answer 1

Since at the equilibrium time derivatives are zero, the issue reduces to solving a system of algebraic equations. Mathematica has many commands for equation solving, among them Solve and Reduce.

Here is a snapshot of solving your particular system:

Answer 2

In Maple, it's pretty straightforward.

eqs:= {mu - (beta*i + mu)*S = 0, beta*S*i - (mu+sigma)*E = 0, sigma*E - (mu+gamma)*i = 0, gamma*i - mu*R = 0}; solve(eqs,{S,E,i,R});

{E = 0, R = 0, S = 1, i = 0}, {E = -mu*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/sigma/(mu+sigma)/beta, R = -gamma*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/(mu^2+mu*gamma+sigma*mu+sigma*gamma), S = (mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/sigma, i = -mu*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/(mu^2+mu*gamma+sigma*mu+sigma*gamma)}