Maximize and FindMaximum drive me crazy

Question

I have to run a numerical study. As part of it I have to solve:

A = 1*10^6 ; b = 1/4; delta = 0; 
K[x_] = 0.9772498680518208` + 
    x (1.186104449210508`*^-10 + 
       x (-1.4649246311182838`*^-19 + 
          x (3.0214683712038235`*^-28 + 
             x (-2.6422357409090615`*^-37 + 
                x (-5.8010340178270696`*^-46 + 
                   x (-2.6339819094686574`*^-55 + 
                    x (-3.7463204038553128`*^-65 + \
(6.191585769652652`*^-77 + 6.879539744126434`*^-87 x) x))))))) // 
   FullSimplify;

F[x_] = Piecewise[{{0, -10^10 < x < -2.2*10^9}, {K[x], .2*10^9 >= 
      x >= -2.5*10^9}, {1, .2*10^9 < x}}];
H1[a_, i_, j_] := 
 NMaximize[{(A - a p) (p - w) +         (1 - 
       F[-(A - a p) (w - i)]) (A - a j)^2/(4 a)           + 
    F[-(A - a p) (w - i)] ( (A - a i)^2/(8 a) - b (A - a p) (p - w)), 
   i <= w <= A/a, 0 < p <= A/a}, {w, p}]
H2[a_, i_, j_] := 
 Maximize[{(A - a p) (p - w) +         (1 - 
       F[-(A - a p) (w - i)]) (A - a j)^2/(4 a)           + 
    F[-(A - a p) (w - i)] ( (A - a i)^2/(8 a) - b (A - a p) (p - w)), 
   i <= w <= A/a, 0 < p <= A/a}, {w, p}]
H3[a_, i_, j_] := 
 FindMaximum[{(A - a p) (p - w) +         (1 - 
       F[-(A - a p) (w - i)]) (A - a j)^2/(4 a)           + 
    F[-(A - a p) (w - i)] ( (A - a i)^2/(8 a) - b (A - a p) (p - w)), 
   i <= w <= A/a, 0 < p <= A/a}, {w, p}]
H1[50, 5000, 6500]
H2[50, 5000, 6500]
H3[50, 5000, 6500]

H1[50, 5000, 5000]
H2[50, 5000, 5000]
H3[50, 5000, 5000]

Sry for the long functions, but perhaps this is important to identify the error.

When I run this code, the first two H1 and H2 find in the first case a smaller maximum than H3, in the second case a higher maximum. I would like to always get the "real", that is, the higher maximum.

Any suggestions?

Best regards Andreas

Answer 1

The answer is "Know thy function."

Because you just have two variables to determine and you know the ranges for those variables, a contour plot will show where there might be some local minima and maxima:

First, modify H3 to allow for starting values:

H3[a_, i_, j_, w0_, p0_] := FindMaximum[{(A - a p) (p - w) + 
  (1 - F[-(A - a p) (w - i)]) (A - a j)^2/(4 a) + 
  F[-(A - a p) (w - i)] ((A - a i)^2/(8 a) - b (A - a p) (p - w)), 
  i <= w <= A/a, 0 < p <= A/a}, {{w, w0}, {p, p0}}]

Create a function to construct a contour plot:

h[a_, i_, j_] := ContourPlot[(A - a p) (p - w) + (1 - F[-(A - a p) (w - i)]) (A - a j)^2/(4 a) + 
  F[-(A - a p) (w - i)] ((A - a i)^2/(8 a) - b (A - a p) (p - w)), 
  {w, i, A/a}, {p, 0, A/a}, PlotPoints -> 100, PlotRange -> All, 
  Contours -> 10^8 Table[c, {c, 10, 55}]]

Now look at the contour plots and associated solutions:

h1 = H1[50, 5000, 6500]
(* {3.45048*10^9, {w -> 8257.58, p -> 12500.}} *)
h2 = H2[50, 5000, 6500]
(* {3.45048*10^9, {w -> 8257.58, p -> 12500.}} *)
h3 = H3[50, 5000, 6500, 6000, 12000]
(* {3.55146*10^9, {w -> 5000., p -> 12500.}} *)
Show[h[50, 5000, 6500],
 ListPlot[{{{w, p} /. h1[[2]]}, {{w, p} /. h2[[2]]}, {{w, p} /. h3[[2]]}},
  PlotStyle -> {{Blue, PointSize[0.03]}, {Green, PointSize[0.02]}, {Red, PointSize[0.01]}},
  PlotLegends -> {"H1", "H2", "H3"}]]

h1 = H1[50, 5000, 5000]
(* {3.85317*10^9, {w -> 8844.11, p -> 12500.}} *)
h2 = H2[50, 5000, 5000]
(* {3.85317*10^9, {w -> 8844.11, p -> 12500.}} *)
h3 = H3[50, 5000, 5000, 1, 1]
(* {3.56361*10^9, {w -> 5000., p -> 12500.}} *)
Show[h[50, 5000, 5000],
 ListPlot[{{{w, p} /. h1[[2]]}, {{w, p} /. h2[[2]]}, {{w, p} /. h3[[2]]}},
  PlotStyle -> {{Blue, PointSize[0.03]}, {Green, PointSize[0.02]}, {Red, PointSize[0.01]}},
  PlotLegends -> {"H1", "H2", "H3"}]]

When using FindMaximum (and FindMinimum good starting values are usually essential.