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

-