Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question

3 Answers 3

up vote 3 down vote accepted

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:

enter image description here

share|improve this answer
    
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)}

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.