RegionPlot won't finish, but Plot3D does

Question

I would like to produce a region plot where some function is less than 1. This runs for hours without finishing (until I lose patience and abort it). However, I can plot the function (using Plot3D), which finishes in a few minutes. I have tried putting in N's, etc., to force it to evaluate numerically to no avail.

In some more detail: the function is naively a function of three variables, $x$, $y$, $z$, but I use a constraint to (numerically) solve for $z$ in terms of $x$ and $y$. This reduces it to a function of two variables, which should be amenable to RegionPlot. I at first thought something was going wrong with the numerical solution for $z$ in terms of $x$ and $y$, but Plot3D is able to evaluate the function on the same $(x,y)$ region just fine.

"Minimal" working example:

Δ[z_, x_, a_] := (z^2*x/2)*(1 - a/(z*x))^2; 
v = (220*10^3)/(3*10^8);
σ[z_, x_] := (4*Pi*z^2)/x^2*(2^10*Pi*z/v)/(3*Exp[4]);  
Γ[z_, x_, y_] := σ[z, x]/y^2*8.440512301344581`*^25* .9*
                              Sqrt[1 - (2*10^-3)^2/(Δ[z, x, 2*10^-3] + x*v^2/2)^2];
zCons[x_, y_] := z /. NSolve[Γ[z, x, y] == 1.2000000000000002`*^19, z, Reals][[1]];
β[z_, x_, y_] := 10^-10/z^4* (y *x);
Plot3D[β[zCons[x, y], x, y], {x, 5, 45}, {y, 150, 350}]
RegionPlot[β[zCons[x, y], x, y] < 1, {x, 5, 45}, {y, 150, 350}]

My apologies that this is still rather complicated, but it's based on the actual function I'm trying to plot. On my computer, it takes about 3 minutes for the Plot3D command to run, but RegionPlot goes on for at least 3 hours. In that time, it only gives the Solve::ratnz error, which is fine; I don't mind it using numerical techniques.

The output from Plot3D is:

which does not show any pathological problems.

Thank you very much for you time and help.

Answer 1

This answer perhaps gets at a workaround, rather than a solution. I find it easier to create a table of discrete data under circumstances where I am plotting number-crunching-heavy functions.

data = Table[{x, y, \[Beta][zCons[x, y], x, y]}, {y, 150, 350, 5}, {x, 5, 45, 5}];

Your region data can be obtained through interpolating the data:

RegionPlot[Interpolation[
    Flatten[data, {2, 1}]][x, y] < 1, {x, 5, 45}, {y, 150, 350}]

Perhaps not as elegant as a one-step solution, but I can make the table and the plot in less time than it took to run Plot3D on the functions themselves. For example:

Module[{}, 
  data = Table[{x, y, \[Beta][zCons[x, y], x, y]}, {y, 150, 350, 
     5}, {x, 5, 45, 5}];
  Plot3D[Interpolation[Flatten[data, {2, 1}]][x, y], {x, 5, 45}, {y, 
    150, 350}]] // AbsoluteTiming
Plot3D[\[Beta][zCons[x, y], x, y], {x, 5, 45}, {y, 150, 
   350}] // AbsoluteTiming

The region plot command required 16 seconds on the same computer.

Answer 2

You can save yourself some time by simplifying the equation NSolve has to deal with. It's the slow piece of your code. NSolve works most efficiently on polynomial equations. It can deal with some other equations, but in this case, it seems to take a route that's not as efficient as what is mathematically obvious.

If we square both sides of the equation, the square root goes away. Next, if we put the terms on one side and combine them, then we just need to solve where the numerator is zero. Squaring both sides introduces an extraneous solution, but the one we want turns out to be the last one (for every x and y).

Here is the code:

neweq = Γ[z, x, y]^2 - (1.2000000000000002`*^19)^2 // Expand // Together // Numerator;
zCons2[x_?NumericQ, y_?NumericQ] := z /. NSolve[neweq == 0, z, Reals][[-1]];

Old:

Plot3D[β[zCons[x, y], x, y], {x, 5, 45}, {y, 150, 350}] // AbsoluteTiming

New:

Plot3D[β[zCons2[x, y], x, y], {x, 5, 45}, {y, 150, 350}] // AbsoluteTiming
RegionPlot[β[zCons2[x, y], x, y] < 1, {x, 5, 45}, {y, 150, 350}] // AbsoluteTiming

Answer 3

First, much thanks to @bobthechemist and @Michael_E2; using their solutions I was able to track down the problem and solve it so that I can use RegionPlot with the original function.

Following @bobthechemist, I evaluated the function on ever finer grids. The times increased but remained reasonable. When the grid got approximately as fine as the mesh used in the Plot3D command, making the Interpolated RegionPlot took about the same amount of time as making the 3D plot.

On the other hand, if I set MaxRecursion -> 0, PlotPoints -> 2 in the RegionPlot (calling the original function), it still never returned, even though it should be a very crude mesh. This suggested to me that for some reason it was having problems evaluating it even at a single point. As @Michael_E2 pointed out, the function is not a nice polynomial, and converting it into one does resolve the issue. However, he seems to suggest that the problem is simply that as written the NSolve is too slow; to convince myself that it wasn't, I let my code run overnight. As expected, it ran until my computer ran out of memory, and looking at the stack showed that it has still not evaluated a single point in the mesh.

As I wrote in my question, I did try using _?NumberQ and N[ ] to try to force it to evaluate the expression inside the RegionPlot command numerically, but a search through the stack confirmed my suspicion (developed after @bobthechemist's grid method worked) that for some reason it was still trying to handle the condition inside the RegionPlot command symbolically, and wasn't getting anywhere.

Since N[ ] didn't work, I decided to give a generous use of Hold/ReleaseHold a try. That seems to do the trick, although I'm not entirely sure why it works where N[ ]'s didn't. Here's code in which the Plot3D and the RegionPlot run in approximately the same amount of time:

Δ[z_, x_, Mγ_] = (z^2*x/2)*(1 - Mγ/(z*x))^2; 
v= (220*10^3)/(3*10^8);
σ[z_, x_] = (4*Pi*z^2)/x^2*(2^10*Pi*z/v)/(3*Exp[4]);  
Γ[z_, x_?NumberQ, y_?NumberQ] = σ[z, x]/
y^2*8.440512301344581`*^25* .9*Sqrt[1 - (2*10^-3)^2/(Δ[z, x, 2*10^-3])^2];
zCons[x_?NumberQ, y_?NumberQ] = Hold[z /. NSolve[Rationalize[Γ[z, x, y] == 1.2000000000000002`*^19, 10^-10], z, Reals][[1]]];
β2[x_?NumberQ, y_?NumberQ] = Hold[10^-10/(ReleaseHold[zCons[x, y]])^4* (y *x)];
RegionPlot[ReleaseHold[β2[x, y]], x, y] < 1, {x, 5, 45}, {y, 150, 350}]