Crthuky

This is a question that follows up from this post:

Select one among multiple solutions that satisfies a certain condition and 3D Plot it with varying simulation values

I was trying to apply the answer to a more complicated case I am working on, and I'm getting an error message.

To start, here is the code of my more complicated function:

eq = 1/(24 (-1 + r)^3 r^4 s) (3 d^6 (-3 + 5 r) (-1 + t) + 6 d^5 (3 - 5 r + c (-3 + r (9 - 4 q + 8 (-1 + q) r))) (-1 + t) + 12 d (-1 + r)^3 r^2 (s + 2 t) - 4 (-1 + r)^3 r^2 s ((-1 + r^2) s (-1 + t) + 3 t) + 12 d^3 (-1 + r) r (2 + c (-2 + r (5 - 2 q + 4 (-1 + q) r)) + s - r (3 + s + r s) + 2 t - 3 r t + (-1 + r + r^2) s t) + 3 d^4 (3 + c^2 (1 + 2 (-1 + q) r) (-3 + r (7 - 2 q + 6 (-1 + q) r)) (-1 + t) - 2 c (-3 + r (9 - 4 q + 8 (-1 + q) r)) (-1 + t) - 3 t + r (3 + 4 r (-5 + 3 r) + 5 t)) + 12 d^2 (-1 + r) r (-(-1 + c) (s (-1 + t) - 2 t) + r^3 (-1 + c (-1 + q) s (-1 + t) - t) + r^2 (2 + (-1 + c (-1 + 2 q)) s (-1 + t) + (2 + 4 c - 4 c q) t) + r (-1 - (1 + c (-2 + q)) s (-1 + t) + (2 - 5 c + 2 c q) t))) == 0

Among multiple solutions of $r$, I would like to select the positive real $r$ and simulate it with the following simulation values, for example, $d=0.8$, $s=2$, $t=0.02$ and $0<c<1$ and $1<q<2$. I would like to 3DPlot $r$ against $c$ and $q$.

My Mathematica code is as follows:

With[d = 0.8, s = 2, t = 0.02,
 sol = Solve[eq, 0 <= c <= 1, 1 <= q <= 2, r, r > 0]]; 
Plot3D[Evaluate[r /. sol], c, 0, 1, q, 1, 2]

And once run, I get a very long warning message including something like:

Warning: r>0 is not a valid domain specification. Assuming it is a variable to eliminate.

Can anyone help?

Mathematica's kernel seemed to slow down to a crawl when constraints on c and q were given along with your equation. However, the constraints are enforced in the s plot domain specification, so I thought it would worth trying to solve the equation without them which succeeded. Rationalization of the parameters suppressed all warnings.

Making the plot was a little slow because there seem some regions with very steep gradients that Plot3D has trouble with.

eq = (1/(24 (-1 + r)^3 r^4 s) (3 d^6 (-3 + 5 r) (-1 + t) + 6 d^5 (3 - 5 r + c (-3 + r (9 - 4 q + 8 (-1 + q) r))) (-1 + t) + 12 d (-1 + r)^3 r^2 (s + 2 t) - 4 (-1 + r)^3 r^2 s ((-1 + r^2) s (-1 + t) + 3 t) + 12 d^3 (-1 + r) r (2 + c (-2 + r (5 - 2 q + 4 (-1 + q) r)) + s - r (3 + s + r s) + 2 t - 3 r t + (-1 + r + r^2) s t) + 3 d^4 (3 + c^2 (1 + 2 (-1 + q) r) (-3 + r (7 - 2 q + 6 (-1 + q) r)) (-1 + t) - 2 c (-3 + r (9 - 4 q + 8 (-1 + q) r)) (-1 + t) - 3 t + r (3 + 4 r (-5 + 3 r) + 5 t)) + 12 d^2 (-1 + r) r (-(-1 + c) (s (-1 + t) - 2 t) + r^3 (-1 + c (-1 + q) s (-1 + t) - t) + r^2 (2 + (-1 + c (-1 + 2 q)) s (-1 + t) + (2 + 4 c - 4 c q) t) + r (-1 - (1 + c (-2 + q)) s (-1 + t) + (2 - 5 c + 2 c q) t))));

Module[eqn, soln,
 eqn = With[d = 4/5, s = 2, t = 2/100, Evaluate @ eq];
 soln = Solve[Evaluate @ eqn == 0, r];
 Plot3D[Evaluate[r /. soln], c, 0, 1, q, 1, 2, PlotPoints -> 50]]

Edit

This addresses a issue raised by the OP in a comment below.

Here is way to get numerical values of r.

soln =
 Module[eqn,
 eqn = With[d = 0.8, s = 2, t = 0.02, Evaluate@eq];
 NSolve[Evaluate @ eqn == 0, r]];
Do[rVal[i][c_, q_] = soln[[i, 1, 2]], i, Length[soln]];

Then to get the value of r for the 1st solution with c = .5 and `q = 1.5:

rVal[1][.5, 1.5]

-0.485047

To get all the value of r for all seven solutions at c = .25 and `q = 1.

Through[Array[rVal, 7][.25, 1.]]

-0.524898, 0.519435, 0.792331, 0.986375, 1.19151, 
 0.0176256 - 0.0442494 I, 0.0176256 + 0.0442494 I

Notice that in this last case there are five real solutions.

2nd Edit

Here is an exploration of the seven root objects. Note that two of plots are are empty because two the root objects are never real in the specified domain. Also, you should be able to see where the icicles come from.

Column[
 Table[
 Plot3D[Evaluate[rVal[i][c, q]], c, 0, 1, q, 1, 2, PlotPoints -> 50], 
 i, 7]]