# Find function tracing the maximum of a 2D function

I'd like to solve the following problem: Given

MI[o_, n_] := Log2[1 + 4 o (Sqrt@n + Sqrt[1 + n])^2]


Find a function f[n_] that tells which o to use to optimise (maximise) MI.

My first pen and paper approach was to calculate the gradient of the function. Then, starting from the origin, the idea is to step "a little bit" into this direction to get to the next, best, point. However, for a closed solution, this means solving an iteration problem. So I went to Mathematica. Here, the approach looks like

Maximize[{MI[o, n], 0 <= Sqrt[o^2 + n^2] <= b, b > 0}, {o, n}]


This is out of Mathematica's (10.2) power. What works at least is

a = 3;
Simplify@Maximize[{MI[o, n], 0 <= Sqrt[o^2 + n^2] <= a}, {o, n}]


If someone has an idea of an approach that leads to a symbolic solution, I'd be happy!

PS: As @march, @rhermans and @Jack LaVigne pointed out, the question only makes sense when using a bound to o and n, which I did not mention (besides implicitly in the last code).

• What makes you believe your expression does have a maximum? – rhermans Nov 12 '15 at 16:53
• You have a constraint on o and n in your Maximize function. Are you maximizing MI subject to this constraint? Because otherwise, as @rhermans notes, it seems like your function grows monotonically with both o and n. – march Nov 12 '15 at 17:07
• If constrained to a disk you could use Maximize[mi[o, n], {o, n} \[Element] Circle[{0, 0}, a], Reals] but apparently there is no simple general form for that. – rhermans Nov 12 '15 at 17:29

If you make a plot of your function over a region you will see that it grows to infinity as o and n increase.

mi[o_, n_] := Log2[1 + 4 o (Sqrt@n + Sqrt[1 + n])^2]

Plot3D[mi[x, y], {x, 0, 5}, {y, 0, 5}, AxesLabel -> Automatic]


So it is certain that you need to place a bound on it.

max = Maximize[{mi[o, n], 0 <= Sqrt[o^2 + n^2] <= 3}, {o, n}]


outputs an expression containing Root. We can get the numerical value from:

N[max]
(* {6.48452, {o -> 2.2254, n -> 2.01186}} *)


and we can plot the function, a cylinder with radius of three and the solution as a red point.

Show[
Plot3D[mi[x, y], {x, 0, 5}, {y, 0, 5}, AxesLabel -> Automatic] /.
EdgeForm[] -> Sequence[EdgeForm[], Opacity[0.4]],
ParametricPlot3D[{3 Cos[theta], 3 Sin[theta], z}, {theta, 0,
90 Degree}, {z, 0, max[[1]]}],
Graphics3D[{
Red,
PointSize -> 0.04,
Point[{o, n, max[[1]]} /. N[max[[2]]]]
}]
]


From this plot it is clear that the solution will lie on the boundary defined by the constraint.

We can use that fact to eliminate one variable.

We will replace o with Sqrt[b^2 - n^2] (tip from Clemens)

When we make a plot of mi[o,Sqrt[b^2 - n^2]] and observe that it has a similar shape (more or less) independent of the parameter b.

Manipulate[
Plot[mi[Sqrt[b^2 - n^2], n], {n, 0, b}, PlotStyle -> {Thick, Red}],
{{b, 2}, 0.5, 10, Appearance -> "Open"}
]


So now the strategy is to locate where the derivative is zero.

sol = Solve[D[mi[Sqrt[b^2 - n^2], n], n] == 0 && b > 0, n, Reals]


results in

{{n -> ConditionalExpression[Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 1],
0 < b < (3 Sqrt[3/2])/4]},
{n -> ConditionalExpression[Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 3],
b > (3 Sqrt[3/2])/4]}}


So for example if you want to find the value at b=2 you can use:

n /. sol /. b -> 2.
(* {Undefined, 1.31088} *)


Analyze the solution

The solution uses the same Root function but the first and third roots depending upon the value of b.

Let's plot n as a function of the radius b.

f1[b_] := Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 1]

f3[b_] := Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 3]

Show[
Plot[f1[b], {b, 0.001, (3 Sqrt[3/2])/4}, PlotStyle -> Blue,
AxesLabel -> {b, n}],
Plot[f3[b], {b, (3 Sqrt[3/2])/4, 3}, PlotStyle -> Red],
PlotRange -> {{0, 3}, {0, 2.1}}
]


n is a nice smooth curve as a function of b.

Function

Last step is to create a function which outputs o, n, mi as a function of b.

parametricFun[b_] := Module[
{
o,
n,
z
},
n = If[b < (3 Sqrt[3/2])/4,
Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 1],
Root[-b^4 + 2 b^2 #1^2 + #1^3 &, 3]
];
o = Sqrt[b^2 - n^2];
z = mi[o, n];
{o, n, z}
]


Now plot the maximum mi as a function of b.

Show[
Plot3D[mi[x, y], {x, 0, 5}, {y, 0, 5}, AxesLabel -> Automatic] /.
EdgeForm[] -> Sequence[EdgeForm[], Opacity[0.4]],

(ParametricPlot3D[{# Cos[theta], # Sin[theta], z},
{theta, 0, 90 Degree}, {z, 0, parametricFun[#][[3]]},
Mesh -> None] /.
EdgeForm[] -> Sequence[EdgeForm[], Opacity[0.6]])
& /@ Range[1, 7],

ParametricPlot3D[parametricFun[b], {b, 0, 7},
PlotStyle -> {Thick, Black}],

Graphics3D[{Red, PointSize -> 0.025,
Point[parametricFun[#] & /@ Range[7]]
}]
]


• note that ToRadicals[ Root[-b^6 + b^8 + (3 b^4 - 4 b^6) #1^2 + (-3 b^2 + 4 b^4) #1^4 + #1^6 &, 1]] give a "closed form" expression. (Not a nice one...) – george2079 Nov 12 '15 at 21:28
• Solving for n instead of o makes the expressions much nicer. It also prevents from the problem encountered at b == 1 in o /. sol /. b -> 1.. – Clemens Nov 13 '15 at 8:08
• @Clemens Thank you for the observation. I have completed the answer by creating the function and incorporated your tip. – Jack LaVigne Nov 13 '15 at 15:07