0
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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"}]]

Contour plot with j equal 6500

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"}]]

Contour plot with j equal 5000

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

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.