Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

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

Search the Archive

Re: Mathematica hangs when solving sys. of equations w/certain parameters

  • To: mathgroup at smc.vnet.net
  • Subject: [mg55619] Re: [mg55586] Mathematica hangs when solving sys. of equations w/certain parameters
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Thu, 31 Mar 2005 01:24:26 -0500 (EST)
  • References: <200503300821.DAA21889@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Sean wrote:
> I am trying to find local maxima of an objective function by solving the FOCs of a function of 4 control variables (pL,pH,phiL,phiH) and 3 parameters (c, pDBS,phiDBS).
> 
> The objective function is 
> 
> f[pL_,phiL_,pH_,phiH_,pDBS_,phiDBS_,c_] = \
> (pL-0.5c*phiL^2)*((pH-pL)/(phiH-phiL) - pL/phiL) + \
> (pH-0.5c*phiH^2)*((pDBS-pH)/(phiDBS-phiH) - (pH-pL)/(phiH-phiL))
> 
> I find the FOCs by differentiating with respect to each of the 4 control variables, and then I solve the system using Solve[] (see details below)
> 
> For certain values of the parameters (c, pDBS, phiDBS), the Solve[] routine works fine. But I can't figure out why Solve[] just seems to hang for other parameter values. I have verified that it's not because I'm dividing by zero or anything, and the points that are causing problems seem to be totally generic.
> 
> For example, everything works fine when (c=0.1, pDBS=5, phiDBS=8) and when (c=0.1, pDBS=7, phiDBS=8), but not when (c=0.1, pDBS=6, phiDBS=8). I happen to know that at these latter parameter values, one of the solutions to the system is (pL=11.25, pH=0, phiL=15, phiH=0), and when I manually plug this solution into the FOCs, everything evaluates to zero as expected. I have also verified that I am not dividing anything by zero at this solution.
> 
> I was wondering whether anybody else has ever experienced problems with Solve[] working for some parameter values but not for others, when there should be no obvious reason why it wouldn't work.
> 
> Thanks for your help! I'm new to this forum, and did a search to make sure that this question hasn't been asked previously, but didn't see anything of this nature. 
> 
> ---------------------------------------------
> Details
> ---------------------------------------------
> focpL[pL_,phiL_,pH_,phiH_,pDBS_,phiDBS_,c_] = D[f[pL,phiL,pH,phiH,pDBS,phiDBS,c],pL]
> focphiL[pL_,phiL_,pH_,phiH_,pDBS_,phiDBS_,c_] = D[f[pL,phiL,pH,phiH,pDBS,phiDBS,c],phiL]
> focpH[pL_,phiL_,pH_,phiH_,pDBS_,phiDBS_,c_] = D[f[pL,phiL,pH,phiH,pDBS,phiDBS,c],pH]
> focphiH[pL_,phiL_,pH_,phiH_,pDBS_,phiDBS_,c_] = D[f[pL,phiL,pH,phiH,pDBS,phiDBS,c],phiH]
> 
> pDBS:= 6.0
> phiDBS:= 8
> c:=0.1
> 
> solution = \ Solve[{focpL[pL,phiL,pH,phiH,pDBS,phiDBS,c]==0, \
> focphiL[pL,phiL,pH,phiH,pDBS,phiDBS,c]==0, \
> focpH[pL,phiL,pH,phiH,pDBS,phiDBS,c]==0, \
> focphiH[pL,phiL,pH,phiH,pDBS,phiDBS,c]==0},{pL,phiL,pH,phiH}]


I think at least a part of the problem is that the first order 
conditions, with denominators cleared, seem to vanish on a two 
dimensional set. I show this as below.

f = (pL-1/2c*phiL^2)*((pH-pL)/(phiH-phiL) - pL/phiL) +
   (pH-1/2c*phiH^2)*((pDBS-pH)/(phiDBS-phiH) - (pH-pL)/(phiH-phiL))

polys = Numerator[Together[Flatten[Outer[D,{f},{pL,phiL,pH,phiH}]]]]

We form a Groebner basis with respect to a term order that will not take 
forever, using as coefficients the set of rational functions in your 
"parameter" variables.

gb = GroebnerBasis[polys, {pL,phiL,pH,phiH},
   MonomialOrder->DegreeReverseLexicographic,
   CoefficientDomain->RationalFunctions];

Now we check to see which variables appear alone in "head terms". If not 
all appear in that way then the solution set is not finite (it has 
dimension at least one, and, typicically, will actually equal the number 
of variables that do not so appear). DistributedTermsList below will 
show explicitly the internal form of these polynomials using exponent 
vectors, and hence it will be easy to see which "pure powers" appear in 
head terms.

In[62]:= InputForm[Map[First,
   First[Internal`DistributedTermsList[gb, {pL,phiL,pH,phiH},
   MonomialOrder->DegreeReverseLexicographic]]]]

Out[62]//InputForm=
{{{0, 2, 0, 0}, c*phiDBS}, {{0, 1, 0, 2}, -(c*phiDBS)},
  {{1, 0, 0, 2}, c*phiDBS}, {{0, 1, 2, 1}, -6*c^2*pDBS*phiDBS^3 +
    3*c^3*phiDBS^5}, {{1, 1, 1, 1}, -24*c*pDBS*phiDBS^2 + 12*c^2*phiDBS^4},
  {{2, 0, 1, 1}, -32*c^2*pDBS*phiDBS^5 + 16*c^3*phiDBS^7},
  {{2, 1, 0, 1}, 12*c^2*pDBS*phiDBS^3 - 6*c^3*phiDBS^5},
  {{3, 0, 0, 1}, 16*c^2*pDBS*phiDBS^5 - 8*c^3*phiDBS^7},
  {{0, 1, 3, 0}, -576*c^2*pDBS*phiDBS^5 + 288*c^3*phiDBS^7},
  {{1, 1, 2, 0}, -576*c^2*pDBS*phiDBS^5 + 288*c^3*phiDBS^7},
  {{2, 0, 2, 0}, 384*c^2*pDBS*phiDBS^6 - 192*c^3*phiDBS^8},
  {{2, 1, 1, 0}, 576*c^2*pDBS*phiDBS^5 - 288*c^3*phiDBS^7},
  {{3, 0, 1, 0}, -64*c^2*pDBS*phiDBS^6 + 32*c^3*phiDBS^8},
  {{3, 1, 0, 0}, -576*c^2*pDBS*phiDBS^5 + 288*c^3*phiDBS^7},
  {{4, 0, 0, 0}, 128*c^2*pDBS*phiDBS^6 - 64*c^3*phiDBS^8}}

Only the first two variables appear alone in head terms. And in fact 
experimentation with the parameter values you give, using NSolve, shows 
that the solution sets are two dimensional.

Recall that I am clearing denominators. If I do not do that (as I'll 
show below) then I do not get into this trouble. But I think the extra 
work needed in forcing those denominators not to vanish may be causing 
the hanging behavior you report.

We'll work with the difficult case.

vars = {pL,phiL,pH,phiH};
rats = Block[{c=1/10,pDBS=6,phiDBS=8},
	Together[Flatten[Outer[D,{f},{pL,phiL,pH,phiH}]]]]
polys1 = Numerator[rats]

Now we get denominators and find all nontrivial factors.

In[81]:= polys2 = Map[First,
   Drop[FactorList[Apply[Times,Denominator[rats]]],1]]
Out[81]= {-8 + phiH, phiH - phiL, phiL}

Next we set up new polynomials that, in effect, force these factors to 
not vanish. We do this by multiplying each by a new "reciprocal" 
variable and setting the product equal to unity.

elims = Array[recip,Length[polys2]];
polys2b = polys2*elims-1

Now we solve this system.

InputForm[NSolve[Join[polys1,polys2b], vars, elims]]

Out[8]//InputForm=
{{pL -> 0., phiL -> -15.000000000000007, pH -> 0., phiH -> 0.},
  {pL -> 2.5500000000003755 - 3.306811152757582*I,
   phiL -> 6.000000000000587 - 4.898979485566834*I,
   pH -> 5.700000000001064 - 9.553009996855195*I,
   phiH -> 12.000000000001176 - 9.797958971133658*I},
  {pL -> 2.5500000000003755 + 3.306811152757582*I,
   phiL -> 6.000000000000587 + 4.898979485566834*I,
   pH -> 5.700000000001064 + 9.553009996855195*I,
   phiH -> 12.000000000001176 + 9.797958971133658*I},
  {pL -> 11.250000000000094, phiL -> 15.000000000000115, pH -> 0.,
   phiH -> 0.}, {pL -> 1.3055555555559843 + 1.614937983749258*I,
   phiL -> 3.3333333333341413 + 2.981423969998686*I,
   pH -> 2.722222222223331 + 4.223683957498015*I,
   phiH -> 6.666666666668278 + 5.962847939997368*I},
  {pL -> 1.3055555555559843 - 1.614937983749258*I,
   phiL -> 3.3333333333341413 - 2.981423969998686*I,
   pH -> 2.722222222223331 - 4.223683957498015*I,
   phiH -> 6.666666666668278 - 5.962847939997368*I}}

Removing the ones with complex values gives

In[9]:= Select[%,FreeQ[#,Complex]&]
Out[9]= {{pL -> 0., phiL -> -15., pH -> 0., phiH -> 0.},
   {pL -> 11.25, phiL -> 15., pH -> 0., phiH -> 0.}}


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: Mathematica hangs when solving sys. of equations w/ certain parameters
  • Next by Date: Re: Much faster ConvexHull implementation
  • Previous by thread: Re: Mathematica hangs when solving sys. of equations w/ certain parameters
  • Next by thread: Re: Mathematica hangs when solving sys. of equations w/ certain parameters