# RegionPlot of the Maximum of a function

For the function fun[a_, b_, x_, y_] = Sin[a x] Cos[b y] (a x + b y);corresponding to $$x$$ and $$y$$, I need a RegionPlot showing the maximum of function $$fun[a, b, x, y]>0.5$$ such that the maximization is carried over $$a,b\in [0,1000]$$.

A failed attempt: RegionPlot[{NMaximize[ fun[a, b, x, y] > 0.5, {a, 0, 1000}, {b, 0, 1000}]}, {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}]

• Is this really that difficult?
– Mike
Jan 30, 2021 at 13:13
• it looks like the region will be the whole x-y region (Rectangle[{0,0},{2Pi, 2Pi}]) if threshold is less than MaxValue[{fun[a,b,x,y], 0<=a<=1000,0<=b<=1000,0<=x<=2Pi,0<=y<=2Pi}, {a,b,x,y}] and empty region otherwise.
– kglr
Jan 31, 2021 at 12:25

Inequalities are boolean expressions; they do not have maximums and minimums.

\$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global*"]

fun[a_, b_, x_, y_] = Sin[a x] Cos[b y] (a x + b y);

NMaximize[{fun[a, b, x, y], 0 <= a <= 1000, 0 <= b <= 1000, 0 <= x <= 2 Pi,
0 <= y <= 2 Pi}, {a, b, x, y}, WorkingPrecision -> 20] // N

(* {2797.59, {a -> 404.339, b -> 595.551, x -> 3.68672, y -> 2.19444}} *)


For nonlinear expressions, NMaximize generally only finds a local maximum.

Using a non-default Method provides a better result:

(sol = {#, NMaximize[{
fun[a, b, x, y], 0 <= a <= 1000, 0 <= b <= 1000, 0 <= x <= 2 Pi,
0 <= y <= 2 Pi},
{a, b, x, y},
WorkingPrecision -> 20,
Method -> #]} & /@
"DifferentialEvolution", "SimulatedAnnealing"} //
SortBy[#, #[[2, 1]] &] &) // N //
Grid[#, Frame -> All] & (param = sol[[4, 2, 2, 1 ;; 2]]) // N

(* {a -> 973.25, b -> 1000.} *)


Due to the complicated nature of the region's structure, this is quite slow.

RegionPlot[(fun[x, y, a, b] /. param) > 1/2,
{x, 0, 2 Pi}, {y, 0, 2 Pi},
PlotPoints -> 75,
WorkingPrecision -> 20] Using a much higher threshold.

RegionPlot[(fun[x, y, a, b] /. param) > 2500,
{x, 0, 2 Pi}, {y, 0, 2 Pi},
PlotPoints -> 75,
WorkingPrecision -> 20] • Write fun= (ax + by) Cos[by] Sin[ax]  with ax==ax, bx==bx, regard the two terms separated, with the help of Maximize[{Sin[ax] Cos[by], 2 Pi 998 < ax < 2 Pi 1000, 2 Pi 998 < by < 2 Pi 1000}, {ax, by}, Reals]  to get the two equal maxima NMaximize[{(ax + by) Cos[by] Sin[ax], 2 Pi 999 < ax < 2 Pi 1000, 2 Pi 999 < by < 2 Pi 1000}, {ax, by}]  {12561.7, {ax -> 6281.61, by -> 6280.04}}  and NMaximize[{(ax + by) Cos[by] Sin[ax], 2 Pi 999 < ax < 2 Pi 1000, 2 Pi 999 < by < 2 Pi 1000}, {ax, by}, Method -> "DifferentialEvolution"]  {12561.7, {ax -> 6278.47, by -> 6283.19}}  Oct 28, 2021 at 9:59
• max is (7997 \[Pi])/2 ` Oct 28, 2021 at 10:02