# Maximizing a function defined using NIntegrate [duplicate]

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?

## marked as duplicate by Michael E2, Oleksandr R., MarcoB, Öskå, Mr.Wizard♦Aug 31 '15 at 22:29

• The problem is that NMaximize evaluates your functions many times, and each evaluation takes a lot of time, because it involves a numerical integration. A possible solution is to pre-calculate fn on a grid of values, construct an InterpolatingFunction from these data, and run the maximization function on that. – yohbs Aug 31 '15 at 16:13
• I need to calculate x and y  till at least 4 decimal places, and then I may have huge grid. Thats why I am looking for any precise in-built function is available in mathematica. Further, I dont know how I can use interpolating function. – Frey Aug 31 '15 at 16:43
• This answer and its links address how to use a numerical function inside another. – Michael E2 Aug 31 '15 at 17:48

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},
EvaluationMonitor :> (AppendTo[ac, {x, y, g[x, y]}])]

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


• @ belisarius, thanks for your great help !!! You code really helps my work, and will be an interesting technique for others. Have a nice day :) – Frey Aug 31 '15 at 22:45

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.

• @ Oliver Jennrich, thanks for your answer, and it helps for me. As @belisarius mentioned, this may be work faster by setting WorkingPrecision. – Frey Aug 31 '15 at 22:47