Re: FindMaximum doesn't converge

*To*: mathgroup at smc.vnet.net*Subject*: [mg112059] Re: FindMaximum doesn't converge*From*: Yaroslav Bulatov <yaroslavvb at gmail.com>*Date*: Sat, 28 Aug 2010 07:01:56 -0400 (EDT)*References*: <i57rmr$1g2$1@smc.vnet.net>

On Aug 27, 1:06 am, Themis Matsoukas <tmatsou... at me.com> wrote: > The local maximum is an artifact in the contour plot. Your function has a flat maximum with respect to j at h=0: > > FullSimplify[obj /. h -> 0] > > -Log[8] > > You can examine your function using > > Manipulate[ > Column[{Dynamic[H], > Plot[obj /. h -> H, {j, -5, 5}, PlotRange -> {-15, 0}]}], {H, -5, > 5}] > > Themis Well, it's not really a maximum at h->0...more like, h=0,j\in {- infinity,infinity} is a level set I got some help from Daniel Lichtblau, and it turns out that you can use Reduce to find the maximum here. There are two zero-gradient points, a local maximum and a saddle-point. Perhaps the presence of saddle-point is what makes things difficult for FindMaximum obj = 2/3 Log[(E^(-(h/Sqrt[3])) + E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j))/(2 E^(-(h/Sqrt[3])) + 2 E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))] + 1/3 Log[(2 E^(-(h/Sqrt[3])) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j))/(2 E^(-(h/Sqrt[3])) + 2 E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))] + 2/3 Log[(2 E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))/(2 E^(-(h/Sqrt[3])) + 2 E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))] + 4/3 Log[(E^(-(h/Sqrt[3])) + E^(h/Sqrt[3]) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))/(2 E^(-(h/Sqrt[3])) + 2 E^(h/Sqrt[3]) + E^(-(h/Sqrt[3]) - Sqrt[2] j) + E^(h/Sqrt[3] - Sqrt[2] j) + E^(-Sqrt[3] h + Sqrt[2] j) + E^(Sqrt[3] h + Sqrt[2] j))]; exprs = Together[{D[obj, h], D[obj, j]}]; sols = List@ToRules[Reduce[exprs == 0 && h >= 0 && j >= 0, {h, j}= ]]; hess = Table[D[obj, {{j, h}, 2}] /. sol, {sol, sols}]; vals = Table[obj /. extr, {extr, sols}] // N Eigenvalues /@ N[hess]