Symbolic solution of a system of 7 nonlinear equations

Question

I've got a system of ordinary differential equations - 7 equations, and ~30 parameters governing their behavior as part of a mathematical model of disease transmission. I'd like to find the steady states for those equations Changing dx/dt = rest of the equation to 0 = equation for each of the equations makes it a straightforward algebra problem. This could be done by hand, but I'm ridiculously bad at that kind of computation.

I've tried using Mathematica, which can handle smaller versions of this problem (see here), but Mathematica is grinding to a halt on this problem. Is there a more efficient/effective way to approach this? A more efficient symbolic math system? Other suggestions?

A few updates (March 21st):

• The goal is indeed to solve them symbolically - the numerical answers are nice but for the moment the end-goal is the symbolic version.
• There is at least one equilibrium. I haven't actually sat down and proved this, but by design it should have at least one trivial one wherein none is infected at the start. There may not be anything besides that, but that would make me as content as anything else.
• Below is the actual set of equations being talked about.

In summary, I'm looking for symbolic expressions for the solutions of a system of 7 quadratic equations in 7 variables.

Answer 1

It looks like the equations you're dealing with are all polynomial after clearing denominators. That's a good thing (transcendental functions are often a little harder to deal with algebraically). However, it's not a guarantee that your equations have a closed-form solution. This is an essential point that many people don't really "get", even if they know it in theory, so it bears restating: there are fairly simple systems of polynomial equations for which there is no way of giving the solutions in terms of ($n$th) roots etc. A famous example (in one variable) is $x^5-x+1=0$. See also this wikipedia page.

Having said that, of course there are also systems of equations that can be solved, and it's worthwhile to check if your system is one of those. And even if your system cannot be solved, it might still be possible to find a form for your system of equations that is simpler, in some sense. For example, find one equation involving only the first variable (even if it cannot be solved algebraically), then a second equation involving only the first and second variable, etc. There are a few competing theories for how to find such "normal forms" of polynomial systems; the most well-known is Groebner basis theory, and a competing one is the theory of regular chains.

In the computer algebra system Maple (full disclosure: I work for them) both of them are implemented. The solve command typically calls the Groebner basis method, I believe, and that quickly grinds to a halt on my laptop. I tried running the regular chains computation and it takes longer than I have patience for but doesn't seem to blow up as badly memory-wise. In case you're interested, the help page for the command I used is here, and here is the code I used:

restart;
sys, vars := {theta*H - rho_p*sigma_p*
       Cp*(Us/N) - rho_d*sigma_d*D*(Us/N)*rho_a*sigma_a*
       Ca*(Us/N) = 0, 
         rho_p*sigma_p*Cp*(Us/N) + rho_d*sigma_d*
       D*(Us/N)*rho_a*sigma_a*Ca*(Us/N) + theta*H = 0, 
         (1/omega)*Ua - alpha*Up - rho_p*psi_p*
       Up*(H/N) - Mu_p*sigma_p*Up*(Cp/N) - 
             Mu_a*sigma_a*Up*(Ca/N) - Theta_p*
       Up + Nu_up*(Theta_*M + Zeta_*D) = 0, 
         alpha*Up - (1/omega)*Ua - rho_a*psi_a*
       Ua*(H/N) - Mu_p*sigma_p*Ua*(Cp/N) - 
             Mu_a*sigma_a*Ua*(Ca/N) - Theta_a*
       Ua + Nu_ua*(Theta_*M + Zeta_*D) = 0, 
         (1/omega)*Ca + Gamma_*Phi_*D + rho_p*psi_p*
       Up*(H/N) + Mu_p*sigma_p*Up*(Cp/N) + 
             Mu_a*sigma_a*Up*(Ca/N) - alpha*Cp - Kappa_*
       Cp - Theta_p*Cp + Nu_cp*(Theta_*M + Zeta_*D) = 0, 
         alpha*Cp + Gamma_*(1 - Phi_)*D + rho_a*psi_a*
       Ua*(H/N) + Mu_p*sigma_p*Ua*(Cp/N) + 
             Mu_a*sigma_a*Ua*(Ca/N) - (1/omega)*
       Ca - Kappa_*Tau_*Ca - Theta_a*Ca + 
             Nu_ca*(Theta_*M + Zeta_*D) = 
     0, Kappa_*Cp + Kappa_*Tau_*Ca - Gamma_*Phi_*
       D - Gamma_*(1 - Phi_)*D - 
             Zeta_*D + Nu_d*(Theta_*M + Zeta_*D) = 0, 
    Us + H + Up + Ua + Cp + Ca + D = 0, 
         Up + Ua + Cp + Ca + D = 0}, {Us, H, Up, Ua, Cp, Ca, D, N, 
    M}:

sys := subs(D = DD, sys):
vars := subs(D = DD, vars):
params := indets(sys, name) minus vars:
ineqs := [theta > 0 , rho_p > 0 , sigma_p > 
       0 , rho_d > 0 , sigma_d > 0 , 
            rho_a > 0 , sigma_a > 0 , 
      omega > 0 , alpha > 0 , psi_p > 0 , Mu_p > 0 , 
            Mu_a > 0 , Theta_p > 0 , Nu_up > 0 , Theta_ > 
       0 , Zeta_ > 0 , psi_a > 0 , 
            Theta_a > 0 , Nu_ua > 0 , Gamma_ > 0 , Phi_ > 
       0 , Kappa_ > 0 , Nu_cp > 0 , 
            Tau_ > 0 , Nu_ca > 0]:
with(RegularChains):
R := PolynomialRing([vars[], params[]]):
sys2 := map(numer, map(lhs - rhs, normal([sys[]]))):
sol := LazyRealTriangularize(sys2,[],map(rhs, ineqs),[],R);

Answer 2

The professional way is to write your equations in a modeling language such as AMPL or GAMS, and solve it with a solver such as IPOPT.

AMPL is a commercial system, but a free student version of AMPL is able to pose problems with up to 300 equations and variables.

If you just want to solve one or a few problems, you can have it solved online freely by using the NEOS server for optimization - simply submit the AMPL description and wait for the answer to be returned to you.

If you need to solve such systems repeatedly as part of a bigger study (e.g., varying the parameters), you should download IPOPT (which is software under a very liberal licence).

Edit: Note that symbolic solutions that are comprehensible are usually restricted to quite small problems - typically the size of a Groebner basis grows explosively with the number of variables or the degree of the polynomials, and the time for processing even more. Thus a waiting time of an hour or more with Mathematica is a sign (though not a proof) that your symbolic solution would be completely incomprehensible. Moreover, evaluating such a long expression is likely to be numerically unstable, so you would need high precision in the evaluation to get meaningful results.

Answer 3

To write out the entire solution is impossible within reason. But here are some equations to reduce the system down a bit:

$U_{S}$ doesn't appear in any equations other than equations 1 and 2. Furthermore, these equations are a dependent set (equation 1 is -1 times equation 2), so equation 1 can be solved for $U_{S}$ in terms of all other variables, and equation 2 can be discarded.

$$ U_S = H \frac{N \theta (\gamma+\zeta)}{C_A K_A + C_p +K_D} $$ where $K_A = \gamma \rho_A \sigma_A+\kappa \rho_D \sigma_D \tau + \rho_A \sigma_A \zeta$ and $K_D = \gamma \rho_p \sigma_p+\kappa \rho_D \sigma_D + \rho_p \sigma_p \zeta$

Equation 7 is linear in all variables, and can be rearranged to solve for $D$:

$$ D = \frac{\kappa ( C_A \tau + C_p ) } {\gamma+\zeta}. $$

The resulting expression suggests that we should try to solve for any remaining variables in terms of $C_{A}$ and $C_{P}$; since we only have 6 independent equations, the best we can do is to reduce the system to one equation in two variables.

Fortunately, adding equations 3 and 5 together yields an equation that's linear in all variables, and can be solved for $U_{A}$ or $U_{P}$. Adding equations 4 and 6 together also yields an equation that's linear in all variables, and can be solved for $U_{A}$ or $U_{P}$ (whichever wasn't solved for when adding equations 3 and 5 together).

At this point, we should have expressions for $U_{A}$ and $U_{P}$ in terms of $H$, $C_{A}$, and $C_{P}$ (because you can eliminate $D$ using the expression above). We've used equations 1, 2, 5, 6, and 7; we'll keep equations 3 and 4 because they're simpler.

We can use either equation 3 or 4 to solve for $H$ in terms of $C_{A}$ and $C_{P}$. Then, making all of the necessary substitutions, the remaining equation should be only in terms of $C_{A}$ and $C_{P}$. The roots of this equation will determine the steady states of the system; it may or may not be possible to find these roots symbolically.

Good luck!

Answer 4

It depends on the structure of your equations.

If you're looking for all steady states of your set of equations, and you can rearrange them as ErikP says into polynomials, you can use methods from real algebraic geometry to calculate all numerical solutions to high precision. Bertini is one such package that I know of, but there are others. I went to a conference at Notre Dame a few years ago where Bertini was used to find steady states of ODEs from chemical kinetics; Bertini was developed at Notre Dame.

Another possibility is to use the methods proposed in "Nonsmooth exclusion test for finding all solutions of nonlinear equations" by M. D. Stuber, V. Kumar, and P. I. Barton, BIT Numerical Mathematics 50(4), 885-917, DOI: DOI: 10.1007/s10543-010-0280-6; these methods do not require the system of equations to be polynomials. Paul Barton is my adviser, and Matt Stuber is a colleague of mine; if you like, I can ask him for the software and send it to you. The paper uses methods from global optimization and interval arithmetic (it cites ArnoldNeumaier's book), as well as Newton's method. The advantage of this method is that it should locate all solutions; the disadvantage is that it's complicated.

In case it isn't clear, ArnoldNeumaier is suggesting that instead of solving $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ directly using something like Newton's method (which will generally work if you give it a good initial guess, sufficiently close to a solution), you solve

$$\min_{\mathbf{x} \in S}\|\mathbf{F}(\mathbf{x})\|,$$

where $S$ is some feasible set defined by constraints on your problem instead of trying to solve. At a very crude level, using a smooth nonlinear programming solver is a lot like using Newton's method, with additional algorithmic sophistication for robustness and performance. IPOPT is a really good piece of software for this purpose; there are a litany of other solvers out there (just look at the available solvers list for GAMS, AMPL, or NEOS. If you choose a method like this one, be aware of a few caveats:

• It will only locate at most one solution at a time. To find additional solutions, you need to add constraints that exclude all of the previous solutions you've found.
• If your optimization problem is nonconvex, to use IPOPT or similar solvers, you will either need a good initial guess, close to a solution of your equations (same basic principle as Newton's method), or a nonconvex optimization solver like BARON, Couenne, Bonmin, etc. You should try every solver you get your hands on, since the performance of each nonconvex nonlinear programming solver is problem-dependent.

Answer 5

I would suggest looking at a homotopy method. While it is not symbolic, it will produce all the solutions of your problem. For an easy library to check out:

http://homepages.math.uic.edu/~jan/PHCpack/phcpack.html