NMaximize and FindMaximum -- Two Different Results

Question

I have the following code in my preamble:

G = 0.01;
β = 1;
ωc = 50;
ϕ = 0;
θ = π/2;
J = 1;

η = Exp[I ϕ] Tan[θ/2];

Clear[ψ]
ψ[α_, χ_] := Exp[I α]*Tan[χ/2];

integralgamma[ω_, τ_] := 4 G ω Exp[-ω/ωc] ((1 - Cos[ω τ])/ω^(2)) Coth[β ω/2]

integraldelta[ω_, τ_] := 4 G ω Exp[-ω/ωc] (Sin[ω τ] - ω τ)/ω^2

mem : δ[τ_] := 
 mem = NIntegrate[
  integraldelta[ω, τ], {ω, 0, Infinity}, 
  MaxRecursion -> 15, PrecisionGoal -> 3]

mem : γ[τ_] := 
 mem = NIntegrate[
  integralgamma[ω, τ], {ω, 0, Infinity}, 
  MaxRecursion -> 15, PrecisionGoal -> 3]

old[τ_] := (Abs[η]/(1 + Abs[η]^2))^(4 J) Sum[
 Abs[η]^(2 m + 2 p) Binomial[2 J, J + m] Binomial[2 J, J + p] Exp[-I δ[τ] (m^2 - 
  p^2)] Exp[-γ[τ] (m - p)^2], 
   {m, -J, J, 1}, {p, -J, J, 1}];


new[α_, χ_, τ_] := (Abs[ψ[α, χ]]/(1 + Abs[ψ[α, χ]]^2))^(2 J)*(Abs[η]/(1 + Abs[η]^2))^(2 J)*
   Sum[Binomial[2 J, J + m] Binomial[2 J, 
    J + p]*(Conjugate[ψ[α, χ]]*η)^(m)*(Conjugate[η]*ψ[α, χ])^(p)*
      Exp[-I δ[τ] (m^2 - p^2)] Exp[-γ[τ] (m - p)^2], 
    {m, -J, J, 1}, {p, -J, J, 1}];

When I run

J1Final = Plot3D[Re[new[α, χ, 1] - old[1]], {α, 0, 2 π}, {χ, 0, π}, 
    PlotPoints -> 20, MaxRecursion -> 0]

I get:

I try and find the global maximum, which I expect to be near 0.2.

When I run,

NMaximize[{f[α, χ, 1], 0 <= α <= 2 π, 0 <= χ <= π}, {α, χ}]

I get

{-0.0336509, {α -> 6.28319, χ -> 3.14159}}

which is clearly not expected. When I run

FindMaximum[f[α, χ, 1], {α, 2}, {χ, 2}]

I get

{0.291708, {α -> 3.14159, χ -> 1.5708}}

as expected. How could it be that NMaximize is finding the global maximum to be less than what FindMaximum is finding it to be?

Please help on this front. I am now confused if my subsequent analysis with the maximized values may be wrong.

Answer 1

For your example, by using higher precision both NMaximize and FindMaximum find the global maximum.

$Version

(*  "10.4.1 for Mac OS X x86 (64-bit) (April 11, 2016)"  *)

G = 1/100;
β = 1;
ωc = 50;
ϕ = 0;
θ = π/2;
J = 1;
η = Exp[I ϕ] Tan[θ/2];

Clear[ψ, δ, γ, old, new, f]
ψ[α_, χ_] = Exp[I α]*Tan[χ/2];

integralgamma[ω_, τ_] = 
  4 G ω Exp[-ω/ωc] ((1 - 
       Cos[ω τ])/ω^(2)) Coth[β ω/2];

integraldelta[ω_, τ_] = 
  4 G ω Exp[-ω/ωc] (Sin[ω τ] - ω τ)/ω^2;

δ[τ_?NumericQ] := δ[τ] = 
  NIntegrate[integraldelta[ω, τ], {ω, 0, Infinity}, 
   MaxRecursion -> 15, WorkingPrecision -> 15]

γ[τ_?NumericQ] := γ[τ] = 
  NIntegrate[integralgamma[ω, τ], {ω, 0, Infinity}, 
   MaxRecursion -> 15, WorkingPrecision -> 15]

old[τ_?NumericQ] := (Abs[η]/(1 + Abs[η]^2))^(4 J) Sum[
    Abs[η]^(2 m + 2 p) Binomial[2 J, J + m] Binomial[2 J, 
      J + p] Exp[-I δ[τ] (m^2 - 
         p^2)] Exp[-γ[τ] (m - p)^2], {m, -J, J, 1}, {p, -J, J, 1}];

new[α_?NumericQ, χ_?NumericQ, τ_?
    NumericQ] := (Abs[ψ[α, χ]]/(1 + 
        Abs[ψ[α, χ]]^2))^(2 J)*(Abs[η]/(1 + 
        Abs[η]^2))^(2 J)*
   Sum[Binomial[2 J, J + m] Binomial[2 J, 
      J + p]*(Conjugate[ψ[α, χ]]*η)^(m)*(Conjugate[η]*ψ[α, χ])^(p)*
     Exp[-I δ[τ] (m^2 - 
         p^2)] Exp[-γ[τ] (m - p)^2], {m, -J, J, 1}, {p, -J, J, 
     1}];

f[α_?NumericQ, χ_?NumericQ, 1] := 
 Re[new[α, χ, 1] - old[1]]

NMaximize[{f[α, χ, 1], 0 <= α <= 2 π, 
   0 <= χ <= π}, {α, χ}, WorkingPrecision -> 15] // N
(*  N used to suppress display of full precision  *)

(*  {0.291708, {α -> 3.14159, χ -> 1.5708}}  *)

FindMaximum[f[α, χ, 1], {α, 2}, {χ, 2}, 
  WorkingPrecision -> 15] // N
(*  N used to suppress display of full precision  *)

(*  {0.291708, {α -> 3.14159, χ -> 1.5708}}  *)

Answer 2

From the documentation:

• If f and cons are linear, NMaximize can always find global maxima, over both real and integer values.
• Otherwise, NMaximize may sometimes find only a local maximum.

Your function is clearly not linear.

It is possible that in the first series of steps NMaximize found a maximum at the edge, and then refined this part. These functions aren't magic, they work according to specific algorithms, this is why human intervention is sometimes needed. In this case, a symbolic approach seems to work best.