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
- References:
- Mathematica hangs when solving sys. of equations w/ certain parameters
- From: Sean <seanchu@stanford.edu>
- Mathematica hangs when solving sys. of equations w/ certain parameters