Find the maximum of a given function

Question

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.

Answer 1

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[3]]]}}, {(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]