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

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.

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:

• This is exactly what I was talking about, as even solving a system of algebraic equations results in errors because of how genuinely terrible I am at solving anything. – Fomite Nov 4 '11 at 1:22
• I understand most of the code above, but could you possibly walk me through the arguments in the Simplify part? Particularly the beta*sigma =/= (gamma+mu) (sigma+mu) bit. – Fomite Nov 4 '11 at 1:44
• @EpiGrad The simplify was there to remove additional requirements identified by Solve, which were $\beta \not= \frac{(\gamma+\mu)(\mu+\sigma)}{\sigma}$. – Sasha Nov 4 '11 at 1:45

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

To find Equilibrium or Steady State of Differential Equations using MATLAB would be,

clc, clear all;
syms S E I R mu beta sigma gamma 'real'
e1 = eval('mu-(beta*I+mu)*S');
e2 = eval('beta*S*I-(mu+sigma)*E');
e3 = eval('sigma*E-(mu+gamma)*I');
e4 = eval('gamma*I-mu*R');
sol=solve(e1,S,e2,E,e3,I,e4,R);
disp([sol.S sol.E sol.I sol.R])

Output;

$$S=(\gamma\mu + \gamma\sigma + \mu\sigma + \mu^2)/(\beta\sigma),\ \\ E=-(\gamma\mu^2 + \mu^2\sigma +\mu^3 - \beta\mu\sigma + \gamma\mu\sigma)/(\beta\sigma^2 + \beta\mu\sigma),\ \\ I=-(\gamma\mu^2 + \mu^2\sigma + \mu^3 - \beta\mu\sigma + \gamma\mu\sigma)/(\beta\mu^2 + \beta\gamma\mu + \beta\gamma\sigma + \beta\mu\sigma),\ \\ R=-(\gamma^2\mu + \sigma\gamma^2 + \gamma*\mu^2 + \sigma\gamma\mu - \beta\sigma\gamma)/(\beta\mu^2 + \beta\gamma\mu + \beta\gamma\sigma + \beta\mu\sigma)$$