# Is it possible to ask Mathematica to give the maximum value of a function over a specific domain?

I have a function $$f(n)$$ over $$n\in\mathbb{N}$$. I want to know its maximum value over a very large domain, say $$n=1,2,3,4,...,10^6$$. Is there a way to ask Mathematica to give me the maximum value of this function?

f [n_]:= Abs[(8 (-1)^(-2 n) ((11 + 5 Sqrt[5]) Cos[(4 n π)/(
1 + Sqrt[5])] + (7 + 3 Sqrt[5]) Cos[(4 n π)/(
3 + Sqrt[5])] + (7 + 3 Sqrt[5]) π Sin[(4 n π)/(
1 + Sqrt[5])] + (3 + Sqrt[5]) π Sin[(4 n π)/(
3 + Sqrt[5])]))/(
n π^3 (-(2 + Sqrt[5]) Sin[2 n π] +
Sin[2 (-2 + Sqrt[5]) n π])^2)]/n;



NMaximize[{f[n], 1 <= n <= 10^6, n ∈ Integers}, n]

(* Out: {63.041, {n -> 2}} *)


You can then get the corresponding symbolic function value using your definition of f[n]:

f[2] // FullSimplify


• It might be worth noting that NMaximize is not guaranteed to return the maximum. Commented Dec 22, 2020 at 6:07
• @charmin How do you know f[2] is not the global maximum? NMaximize tries to find the global maximum and it might be succeeding (or might not). Commented Dec 22, 2020 at 15:38
• @charmin Yes, I think it matters whether the first local maximum is in fact the global maximum. -- So you need to ask another question? Commented Dec 22, 2020 at 15:46
• @charmin Unfortunately, as far as I know, there is no automatic and general approach to non-linear optimization that guarantees to find the global maximum. Commented Dec 22, 2020 at 16:17
• @MarcoB — well, there is such an approach for integers over a finite domain: exhaustive search. See Michael E2's answer. Commented Dec 22, 2020 at 16:28

Exhaustive search is a reasonable method on a small set:

vals = f[N[Range[10^6], 16]]; // AbsoluteTiming
Position[vals, #] -> # &@Max[vals]
(*
{60.0485, Null}
{{2}} -> 63.04102181896
*)


Confirm f[2] is the max (with arbitrary-precision, we get evidence that round-off error has not messed up the search):

Nearest[vals, Max[vals], 2]
(*  {63.04102181896, 56.03221153}  *)

MinMax[Precision /@ vals]
(* {3.15971, 14.2151}  <-- results have > 3 digits precision *)


If you're sure round-off is no problem, machine-precision is faster:

vals = f[N[Range[10^6]]]; // AbsoluteTiming
Position[vals, #] -> # &@Max[vals]
(*
{1.60528, Null}
{{2}} -> 63.041
*)


Also, using Simplify[f[n], n \[Element] Integers] is 30-35% faster in this case, but that's because f[n] happens to contain terms like Sin[2 n \[Pi]].

• @charmin Don't change the question after it has been answered! If the approach Michael and I have proposed does not work for that new function, then ask a different question. Commented Dec 22, 2020 at 16:16
• @MarcoB I did not mean to change the question, since I could not write my function in the comment, I wrote in the question, then I wanted to remove it. Commented Dec 22, 2020 at 16:32
• Note that depending on the function, you may need to apply Chop or Re to eliminate small imaginary parts that arise in the calculation ((-1)^100000. is not quite real according to Mathematica.) Commented Dec 22, 2020 at 16:38
• @MichaelSeifert Thanks, I forgot that. With arbitrary-precision inputs, the round-off error bounds on the imaginary parts are normally greater than the round-off error bounds, and one can be more rigorous than Chop or Re: E.g. vals /. z_ /; Im@z == 0 :> Re[z]. OTOH, an imaginary part might indicate that the input was not in the domain: E.g., Log[x] - Log[x + 1] /. x -> -2.16, which Chop or Re or my fix would cover up. Perhaps a check like VectorQ[vals, DeveloperRealQ] should precede finding Max[vals]. Commented Dec 22, 2020 at 17:05