WOLFRAM

Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Simultaneous equation problem - too many variables?

Richard Anderson wrote:
> 
> Greetings everyone,
> 
> I have a problem in population genetics that yields the following
> equations:
> 
> Out[84]//InputForm=
> us = (u*(1 + 2*v + h*t*(-1 + \[Theta])) +
>     v*(1 + h*t*(-1 + \[Theta]) - t*\[Theta]))/
>   (2 + u + v + t*(-2 + h*(u + v)*(-1 + \[Theta]) -
>        v*\[Theta]))
> 
> Out[82]//InputForm=
> vs = (u*(1 + 2*v + k*s*(-1 + \[Theta])) +
>     v*(1 + k*s*(-1 + \[Theta]) - s*\[Theta]))/
>   (2 + u + v + s*(-2 + k*(u + v)*(-1 + \[Theta]) -
>        v*\[Theta]))
> 
> u and v are variables. k,h,s,t, and \[Theta] are constants (0 <= k <= 1,
> 0 <= h <= 1, -1 < s < 1, -1 < t < 1). I've tried using
> SetAttributes[s,Constant] and so on, but that doesn't help with the
> following problem.
> 
> I am interested in the solutions of the simultaneous equations {us == u,
> vs == v}. (Please note that us and vs are separate variables, not u*s
> and v*s). A previous question to MathGroup yielded the following idea
> which works for simpler cases of the same kind of problem (the method
> puts (us-u) and (vs -v) over a common denominator so that the problem
> becomes one of solving for the numerator = 0).
> 
> equations ={us==u,vs==v};
> polys=Map[#[[1]]-#[[2]]&,equations]; newpolys=
> Map[Numerator[Together[#]]&,polys]; Reduce[newpolys==0,{u,v}];
> 
> However, in this case, Mathematica(v.3.0.1 running on a Pentium 200 MHz
> and Win 95), runs out of memory, presumably because there are too many
> variables. (As a concession, I can justifiably allow h=k, but if
> possible would like to avoid it. Even with h=k, Mathematica fails to
> solve) I've tried a number of substitutions - e.g. A = 1 + k s(-1 +
> \[Theta]), B = A - t \[Theta] (and similarly in the equation for vs),
> to no avail.
> 
> Any ideas? If anyone needs more information, please e-mail me at the
> address below.
> 
> Thank you.
> 
> Cheers,
> 
> Richard
> 
> ********************************************************************************
> Richard Anderson
> e-mail : richardj.anderson@stonebow.otago.ac.nz


Perhaps because of the intermediate swell, Solve cannot get the system
triangualarized directly. This is a deficiency that I hope to at least
in part address in a future release. But one can use GroebnerBasis,
with options set as per guru specifications, to get the polynomials
into better form.

In[15]:= exprs = {(u*(1 + 2*v + h*t*(-1 + theta)) +
        v*(1 + h*t*(-1 + theta) - t*theta))/
  (2 + u + v + t*(-2 + h*(u + v)*(-1 + theta) -
        v*theta)) - u,
  (u*(1 + 2*v + k*s*(-1 + theta)) +
        v*(1 + k*s*(-1 + theta) - s*theta))/
  (2 + u + v + s*(-2 + k*(u + v)*(-1 + theta) -
        v*theta)) - v};

In[16]:= vars = {u, v};

Now treat all the parameters as part of the field of coefficients. This
is likely to get us a smaller resulting set of polynomials.

In[17]:= Timing[gb = GroebnerBasis[Numerator[Together[exprs]], vars,
        CoefficientDomain->RationalFunctions];] Out[17]= {5.8 Second,
Null}

Let's see how it looks.

In[18]:= Length[gb]
Out[18]= 2

In[19]:= LeafCount[gb]
Out[19]= 3017

In[20]:= Map[Exponent[#,u]&, gb]
Out[20]= {0, 1}

In[21]:= Map[Exponent[#,v]&, gb]
Out[21]= {4, 3}

So we have a univariate quartic polynomial in v, and a linear polynomial
in u with coefficients that are at worst cubics in v. This means in
principle we can solve for v using radicals or Root objects, then back
substitute to solve for u. Now with polynomials of this size you do not
want to attempt radical solutions. But even if I inhibit Roots from
using the complicated radical formulas, it still appears to hang.

SetOptions[Roots, Cubics->False, Quartics->False] Timing[sol =
Solve[gb==0, vars, Sort->False];] (* I abortted this after several
minutes *)

Here are a few other things you can do with gb: (i) Substitute values
for the parameters. (ii )Substitute intervals for the paramters (or a
mix of intervals and explicit values).
(iii) For any fixed theta, minimize some function of u and v, e.g.
u^2+v^2, subject to the parameter constraints on {h,k,s,t}. (iv) Wait a
few years until you have a machine 10x faster and with a 100x the
memory, call Solve (go for the radical solutions), print the result,
and wall-paper a barn with it. Possibly there will be paper left over
so plan accordingly. (Certainly not speaking for my employer on this
one, we are more forest-friendly than this.)

Note that using Simplify is not out of the question, despite the size of
the expressions under scrutiny. For example, if we work with the
special case where h==k, we get a much smaller expression.

In[10]:= gb2 = gb /. h->k;

In[11]:= LeafCount[gb2]
Out[11]= 1927

In[12]:= gb3 = Simplify[gb2];
                                         In[13]:=  LeafCount[gb3]
Out[13]= 420

Even if we do not use the special case we get considerable contraction.
                                                            In[14]:= gb4
= Simplify[gb];                                           

In[15]:= LeafCount[gb4]
Out[15]= 977                                                           

Alternatively we can group the coefficients in {u,v}, then Simplify them
separately.

In[16]:= colgb = Map[Collect[#,{u,v},Simplify]&, gb];

In[17]:= LeafCount[colgb]
Out[17]= 795

Quite likely one of these shorter forms will be more useful in further
manipulations than the raw gb.


Daniel Lichtblau
Wolfram Research
  • Prev by Date: Re: Question about Mathematica
  • Next by Date: RE: Introducing: Conix 3D Explorer!
  • Prev by thread: Simultaneous equation problem - too many variables?
  • Next by thread: Re: Simultaneous equation problem - too many variables?