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]