# Solving an equation with constraints

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 and $$1. 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?

• You should write: Solve[{equations, 0 <= c <= 1, 1 <= q <= 2, r > 0}, r]. As the error suggests, r > 0 does not belong where you have it in your code. – MarcoB Apr 2 at 22:44
• d=0.8;s=2;t=0.02;NMinimize[{(expr)^2,0<c<1&&1<q<2},{r,c,q}] rapidly finds a solution near {r->0.6457, c->0.1130, q->1.5452} Hopefully that will help you, but it isn't clear to me how you intend to Plot3D that. – Bill Apr 2 at 22:57
• @Bill: I'm trying to do simulation to see how the solution for r changes for different values of c and q. Would this be possible? In your code, what does (expr)^2 mean? – ppp Apr 3 at 0:24
• @MarcoB: Would it be possible to simulate (3D Plot) to see how the solution for r changes for different values of c and q? By the way, I'm trying your suggestion, and Mathematica has been running for quite a while. Did you get the answer quickly? – ppp Apr 3 at 0:27
• (expr)^2 was the square of your really long expression but with the trailing ==0 removed. So you wanted to find where your expression equaled zero and this found where the square of your expression was almost exactly zero. I am not sure there are still solutions for expr==0 if you make tiny changes in c or q or both. There may be some solution in r somewhere else, but not close by. I don't know. I'm guessing that is what you are supposed to learn from your exercise. – Bill Apr 3 at 1:32

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 to be 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 an 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 values 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}]]


• Thanks so much! Can you please give me some comments on how to interpret the graphical result? It looks very weird... Why it has layers? What are those icicle-like things? Why change in the color, say from green to purple? Etc... From the graph, what is the approximate value of r for, say, c=0.2 & q=1.2? Isn't it impossible to tell since there are three corresponding values of r? Green, blue, and orange? – ppp Apr 3 at 1:12
• @ppp. Solve finds seven root functions for your equation, so at any patch of the plot domain, up to seven surfaces could show up. But it seems that only three of the root functions actually evaluate to real values over any significant part of the domain. The icicle-like things are areas where the solver found small patches, evidently with very steep gradients, of a real surface for some of the other possible surfaces.There is also a smallish patch of a red surface which can only be seen when the plot is rotated to different view. – m_goldberg Apr 3 at 1:31
• @ppp. You can can explore each of the seven root objects returned from Solve individually to find the one that best represents a solution to your problem. Further, you can get a value of r for any pair of c and q from the individual root objects. – m_goldberg Apr 3 at 1:40
• Do you know how to 'explore each of the root objects returned from Solve individually to find the one that best represents a solution to my problem'? By the way, I have restricted the roots to 0<=r<=1 by adding PlotRange->{0,1}`. It would be desirable to get a unique positive root lying between 0 and 1. – ppp Apr 3 at 1:52
• One more question: I'm sorry but I'm not quite understanding your comments about the icicle-like things. May I ask your explanation once more? I really appreciate! – ppp Apr 3 at 2:05