# Find the maximum of a given function

I calculate a target function f about two variates a and b.

f[a, b] = (1 + k + 2 Sqrt[k] Cos[m]) Sin[a]^2 Sin[2 b]^2 - 4 Sqrt[k]Sin[a/2]^2 Sin[a] Sin[m] Sin[2 b] Sin[4 b] + (1 + k - 2 Sqrt[k] Cos[m]) Sin[a/2]^4 Sin[4 b]^2;


in which 1<k and 0<m<Pi.

And I want to find the maximum while 0<a<Pi and 0<b<Pi/2. The k and m are given and they are real numbers. So I should find the a and b to make the function max. The answer {a,b} should be an expression of k and m. I think this is a problem about finding maximum for a bivar function, so I have tried to calculate the patial derivative and let it be 0 to solve.

Solve[{D[f[a, b], a] == 0, D[f[a, b], b] == 0}, {a, b}, Assumptions -> 0 < a < Pi]


But maybe it is too complex, Mathematica does not give the answer.

Also I use the Reduce function. It can really solve it but it can not be added any assumption so the result is complex and I can not extract the one or two which satisfies the ranges. I also try to use the 2D plot which is about f and b. And observe the influence of the change of a,k,m. But I can not get the analytic solution.

Manipulate[Plot[(1 + k + 2 Sqrt[k] Cos[m]) Sin[a]^2 Sin[2 b]^2 - 4 Sqrt[k]Sin[a/2]^2 Sin[a] Sin[m] Sin[2 b] Sin[4 b] + (1 + k - 2 Sqrt[k] Cos[m]) Sin[a/2]^4 Sin[4 b]^2, {b, 0, Pi/2}], {a, 0, Pi}, {k, 1, 10}, {m, 0, Pi}]


I want to ask if there be any method to get the analytic solution.

• Its numerical solution by f[k_, m_] := NMaximize[{(1 + k + 2 Sqrt[k] Cos[m]) Sin[a]^2 Sin[2 b]^2 - 4 Sqrt[k] Sin[a/2]^2 Sin[a] Sin[m] Sin[2 b] Sin[ 4 b] + (1 + k - 2 Sqrt[k] Cos[m]) Sin[a/2]^4 Sin[4 b]^2, a > 0 && a < Pi && b > 0 && b < Pi && k > 1 && m > 0 && m < Pi}, {a, b}] and the plots Plot3D[f[k, m][[2, 1, 2]], {k, 1, 2}, {m, 0, Pi}] and Plot3D[f[k, m][[2, 2, 2]], {k, 1, 2}, {m, 0, Pi}] and Plot[f[4, m][[2, 2, 2]], {m, 0, Pi}] suggest that optimal values of a and b do not depend on k and these are piece-wise functions of m. Jun 2, 2022 at 6:10
• I find exact solutions art for art's sake. Jun 2, 2022 at 6:10
• But actually is there any analytic solution for a and b about m, rather than the numerical solution. Can we get it? Jun 2, 2022 at 7:18

First numerical analysis give a strong hint, that maximum does not depend on m. Test it with different k and m. Leave it to you to proof it exactly.

f[a_, b_, k_, m_] = (1 + k + 2 Sqrt[k] Cos[m]) Sin[a]^2 Sin[2 b]^2 -
4 Sqrt[k] Sin[a/2]^2 Sin[a] Sin[m] Sin[2 b] Sin[
4 b] + (1 + k - 2 Sqrt[k] Cos[m]) Sin[a/2]^4 Sin[4 b]^2;

nmax[k_?NumericQ, m_?NumericQ] :=
NMaximize[{f[a, b, k, m], {0 <= a <= Pi, 0 <= b <= Pi/2}}, {a, b}]

Plot3D[First@nmax[k, m], {k, 1, 10}, {m, 0, Pi}, PlotRange -> All,
PlotPoints -> 11, MaxRecursion -> 1]


Maximize doesn't do the job for general k, but for given k and arbitrary m. Compare it with nmax[k,m]. Use the trick to substitute an exotic numerical constant for k, evaluate and resubstitute. (Maximize knows the rules to find maximum for any constants as k, but checking routines are not mighty enough to check validity for general k.)

max = Maximize[{f[a, b, k, 0],
0 <= a <= Pi && 0 <= b <= Pi/2 && 1 < k}, {a, b}]

{max = Maximize[{f[a, b, 4, 0], 0 <= a <= Pi && 0 <= b <= Pi/2}, {a,
b}] // FullSimplify, max // N, nmax[4, 0]}

(*   {{9, {a -> \[Pi]/2, b -> \[Pi]/4}}, {9., {a -> 1.5708,
b -> 0.785398}}, {9., {a -> 1.5708, b -> 0.785398}}}   *)

{max = Maximize[{f[a, b, 10 EulerGamma, 0],
0 <= a <= Pi && 0 <= b <= Pi/2}, {a, b}] /. EulerGamma -> k/10 //
FullSimplify, max // N, nmax[10 EulerGamma, 0]}

(*   {{(1 + Sqrt[k])^2, {a -> \[Pi]/2,
b -> \[Pi]/4}}, {(1.\[VeryThinSpace]+ Sqrt[k])^2, {a -> 1.5708,
b -> 0.785398}}, {11.5772, {a -> 1.5708, b -> 0.785398}}}   *)


Different m gives same maximum, but different positions.

{max = Maximize[{f[a, b, 10 EulerGamma, Pi/2],
0 <= a <= Pi && 0 <= b <= Pi/2}, {a, b}] /. EulerGamma -> k/10 //
FullSimplify, max // N, nmax[10 EulerGamma, Pi/2]}

(*   {{(1 + Sqrt[k])^2, {a -> (2 \[Pi])/3,
b -> ArcTan[Sqrt[2 + Sqrt]]}}, {(1.\[VeryThinSpace]+ Sqrt[
k])^2, {a -> 2.0944, b -> 1.09314}}, {11.5772, {a -> 2.0944,
b -> 1.09314}}}   *)


So the function of k for maximum is (1 + Sqrt[k])^2 . Difference to numerical result is as low as nearly workingPrecision.

Plot3D[(1 + Sqrt[k])^2 - First@nmax[k, m], {k, 1, 10}, {m, 0, Pi},
PlotRange -> All]

• +1 Interesting that only m determines the position and only k determines the maximum.
– JimB
Jun 2, 2022 at 18:13