Re: Mathematica hangs when solving sys. of equations w/certain parameters
- To: mathgroup at smc.vnet.net
- Subject: [mg55627] Re: [mg55586] Mathematica hangs when solving sys. of equations w/certain parameters
- From: Chenghuan Sean Chu <cschu at stanford.edu>
- Date: Thu, 31 Mar 2005 01:25:10 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Mr Lichtblau:
Thank you very much for the detailed response. I really appreciate your
help. I'm on a tight deadline and am really swamped with other work at
the moment, but I'll read what you wrote and get back to you.
Regards,
Sean
On Wed, 30 Mar 2005, Daniel Lichtblau wrote:
> 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
>
>
________________________
C. Sean Chu
Stanford University
Landau Economics Building
579 Serra Mall
Stanford, CA 94305
650-248-9648