# Maximize and FindMaximum drive me crazy

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

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[]}, {{w, p} /. h2[]}, {{w, p} /. h3[]}},
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[]}, {{w, p} /. h2[]}, {{w, p} /. h3[]}},
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.