Maximizing a function defined using NIntegrate

Question

I have two functions fn[x,y,z] which is in integral form and I have another function gn[x,y] which is also in integral form included with fn[x,y,z] as well.

m1[x_] := (2 (x - 1))/x;
m2[x_, y_] := (x (2 - x y))/(2 (x - 1) y);

fn[x_, y_, z_]:=
     NIntegrate[E^(-(z/(t m1[x])) - t/m2[x, y])/(t m1[x] m2[x, y]), {t, 0, 1}] + 
     E^(-(1/m2[x, y]))/m1[x] - (E^(-(z/m1[x]) - 1/m2[x, y]) (-1 + E^(z/m1[x])))/m1[x];

gn[x_, y_] := y/2 NIntegrate[Log[1 + z] fn[x, y, z], {z, 0,\[Infinity]}];

I want to find x and y which maximize gn[x,y]. I tried following to calculate x and y but they did not give solutions as they are keep on running.

x /. {#2 & @@ NMaximize[{gn[x, y], 2 > x > 1,1 > y > 0}, {x, y}]}[[1, 1]];
y/. {#2 & @@ NMaximize[{gn[x, y], 2 > x > 1,1 > y > 0}, {x, y}]}[[1, 2]];

Can someone help me to find x and y which maximize gn[x,y] using an alternative method?

Answer 1

Your problem is partly that NIntegrate calls a function that at the given point in time is not full numerical (the first expression in fn). A way around that is to define the functions in such a way that hey will only evaluate for purely numerical values:

Clear[m1, m2, fn, gn];
m1[x_?NumericQ] := (2 (x - 1))/x;
m2[x_?NumericQ, y_?NumericQ] := (x (2 - x y))/(2 (x - 1) y);

fn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := 
NIntegrate[
  E^(-(z/(t m1[x])) - t/m2[x, y])/(t m1[x] m2[x, y]), {t, 0, 1}] + 
  E^(-(1/m2[x, y]))/
  m1[x] - (E^(-(z/m1[x]) - 1/m2[x, y]) (-1 + E^(z/m1[x])))/m1[x];

gn[x_?NumericQ, y_?NumericQ] := 
  y/2 NIntegrate[Log[1 + z] fn[x, y, z], {z, 0, \[Infinity]}];

Still, the numerics are relatively heavy, so some simplification/approximation might be required.

Answer 2

Setting WorkingPrecision -> 5 in the gn integral gives you a reasonable convergence time. At the expense of some more computation time you can check that methods DifferentialEvolution and SimulatedAnnealing both return the sme result given here up to four decimal places.

m1[x_] := 2 (x - 1)/x;
m2[x_, y_] := x (2 - x y)/(2 (x - 1) y);

fn[x_?NumericQ, y_?NumericQ, z_?NumericQ] := 
  fn[x, y, z] = 
   NIntegrate[ E^(-(z/(t m1[x])) - t/m2[x, y])/(t m1[x] m2[x, y]), {t, 0, 1}] + 
    E^(-(1/m2[x, y]))/ m1[x] - (E^(-(z/m1[x]) - 1/m2[x, y]) (-1 + E^(z/m1[x])))/m1[x];

g[x_?NumericQ, y_?NumericQ] := 
  g[x, y] = y/2 NIntegrate[Log[1 + z] fn[x, y, z], {z, 0, ∞},   WorkingPrecision -> 5];

(*The following plot is expensive to calculate, but you may omit it 
because it's done only for displaying *)
plot = Plot3D[g[x, y], {x, 1.1, 2}, {y, .5, 1},   PlotStyle -> {Gray, Opacity[.5]}];

ac = {};
Dynamic@Show[ plot,
             Graphics3D[{PointSize[Large], 
                        Table[{ColorData["SolarColors"][i/Length@ac], Opacity[.5], 
                                Point@ac[[i]]}, {i, Length@ac}]}]
  ]
NMaximize[{g[x, y], 2 > x > 1.1, 1 > y > 0.5}, {x, y}, 
 Method -> {"NelderMead"}, 
 EvaluationMonitor :> (AppendTo[ac, {x, y, g[x, y]}])]

(* {0.13979, {x -> 1.61673, y -> 0.819674}}*)