# Why does not ContourPlot consider the condition given in RegionFunction in the given code?

I have two functions $$f(x,y)$$ and $$g(x,y)$$. I am using  RegionFunction in  ContourPlot to see those values of $$(x,y)$$ for which $$f=g$$ assuming $$-1. I use this code

f := Cos[x y] + (2 Sin[x y])/x + (-(1/2) + Cos[6 x]) Csc[6 x] Sin[x y];
g := Cos[x y] + Sin[x y]/( 2 x) + (-(1/2) + Cos[6 x]) Csc[6 x] Sin[x y] +  Sqrt[-1 + (Cos[x y] + ((1 + 2 x Cot[6 x] - x Csc[6 x]) Sin[x y])/( 2 x))^2];

p1 = ContourPlot[ f == g , {y, 0, 2}, {x, 0, 3},  ContourStyle -> Directive[Black, Thickness[0.006]], PlotPoints -> 50, RegionFunction -> Function[{x, y},  -1 < g < 1   ]];


Then, I check to see if the result of the given plot is correct (i.e. they satisfy the condition $$-1)

p2 = RegionPlot[  -1 < g < 1  , {y, 0, 2}, {x, 0, 3}, PlotStyle -> Green, BoundaryStyle -> Directive[Opacity[0]],   PlotPoints -> 80, FrameLabel -> {"y", "x"}];

Show[{p2, p1}, PlotRange -> All]



and I get this

which shows that  ContourPlot does not consider the condition in  RegionFunction exactly (for some parts). How can I solve this issue? I mean how can I modify the given code to get only those lines which are located in the green area of the given plot?

Reasonable results can be obtained with a few changes to code. First, we can define $$f$$ and $$g$$ like this

ClearAll[f, g]
f[x_, y_] :=
Cos[x y] + (2 Sin[x y])/x + (-(1/2) + Cos[6 x]) Csc[6 x] Sin[x y]
g[x_, y_] :=
Cos[x y] + Sin[x y]/(2 x) + (-(1/2) + Cos[6 x]) Csc[6 x] Sin[x y] +
Sqrt[-1 + (Cos[
x y] + ((1 + 2 x Cot[6 x] - x Csc[6 x]) Sin[x y])/(2 x))^2]


We notice the $$\csc(6x)$$ will produce singularities at $$x=n\pi/6$$, so we want to look at the region $$0 first. Later we will expand this region. Let's start with $$n=1$$, a parameter that is only used in the plot limits.

For the plots we will make $$x$$ horizontal and $$y$$ vertical to avoid confusion. We can change them after we get reasonable results. We must also change the RegionFunction a little. It takes 3 parameters, and they must be in the correct order. The original region function produced a few error messages, but this version seems to be message-free.

n = 1; ymax=6;

p1 = ContourPlot[f[x, y] == g[x, y], {x, 0, n*\[Pi]/6}, {y, 0, ymax},
ContourStyle -> Directive[Black, Thickness[0.006]],
PlotPoints -> 50,
RegionFunction ->
Function[{x, y, z},
x != 0 && Im[g[x, y]] == 0 && -1 < g[x, y] < 1]];

p2 = RegionPlot[
Im[g[x, y]] == 0 && -1 < g[x, y] < 1, {x, 0, n*\[Pi]/6}, {y, 0, ymax},
PlotStyle -> Green, BoundaryStyle -> Directive[Opacity[0]],
PlotPoints -> 80, FrameLabel -> {"x", "y"}];

Show[{p2, p1}, PlotRange -> All, ImageSize -> Small]


Now the value of $$n$$ can be increased and other adjustments can be made.

• Thanks. What is the role of the $z$ variable in your code? @LouisB Commented Feb 14, 2022 at 13:38
• @math2021 ContourPlot passes three values to the RegionFunction, so I put in 3 variable names as arguments. It's an old habit I picked up to check the docs and see what the arguments are. Maybe it's not required, but I would still put it in, if only remind myself that it is available., even if it is not actually used, as in this case. I believe it would have the (Boolean) value of f[x,y] == g[x,y]. Commented Feb 14, 2022 at 14:04