3D Plot: Number of Roots in x of a polynomial in x, a, b and c

I have a polynomial in four variables x,a,b and c. The number of roots of the polynomial in x depends of the choice of a, b and c. I would like to have a 3D-Plot with a, b and c on the axes, while the number of roots >0 at a point (a,b,c) is symbolized by different colours.

The most important to me is to see where exactly the transitions are.

I think I need CountRoots[Polynom,{x,0,Infinity}], but as I am new to Mathematica I can't figure out how to do this. Thanks a lot for your help!

The function is

-a b (1 +
c) (a^8 (1 + c)^4 (6 - 56 c + 68 c^2 - 56 c^3 + 6 c^4 +
3 b (-1 + c)^4 (1 + c)) +
a^7 b (1 + c)^4 (9 b (-1 + c)^4 (1 + c) - 256 c (1 - c + c^2)) +
3 a^6 b^3 (-1 + c)^4 (1 + c)^5 -
8 a^6 b^2 (1 + c)^4 (15 + 52 c - 22 c^2 + 52 c^3 + 15 c^4) -
15 a^5 b^4 (-1 + c)^4 (1 + c)^5 -
128 a^5 b^3 (1 + c)^4 (3 + 2 c + 4 c^2 + 2 c^3 + 3 c^4) +
3 a^3 b^6 (1 + c) (-1 + c^2)^4 -
128 a^3 b^5 (1 + c)^4 (3 + 2 c + 4 c^2 + 2 c^3 + 3 c^4) -
15 a^4 b^5 (-1 + c)^4 (1 + c)^5 -
20 a^4 b^4 (1 + c)^4 (27 + 4 c + 50 c^2 + 4 c^3 + 27 c^4) +
3 a b^8 (-1 + c)^4 (1 + c)^5 -
256 a b^7 c (1 + c)^4 (1 - c + c^2) +
9 a^2 b^7 (-1 + c)^4 (1 + c)^5 -
8 a^2 b^6 (1 + c)^4 (15 + 52 c - 22 c^2 + 52 c^3 + 15 c^4) +
2 b^8 (1 + c)^4 (3 - 28 c + 34 c^2 - 28 c^3 + 3 c^4)) -
a b (1 + c) (21 a^6 b (-1 + c)^4 (1 + c)^3 -
30 a^6 (1 + c)^2 (5 - 12 c + 30 c^2 - 12 c^3 + 5 c^4) +
21 a^5 b^2 (-1 + c)^4 (1 + c)^3 -
60 a^5 b (1 + c)^2 (7 - 4 c + 42 c^2 - 4 c^3 + 7 c^4) -
42 a^3 b^4 (-1 + c)^4 (1 + c)^3 -
120 a^3 b^3 (1 + c) (1 + 37 c + 42 c^2 + 42 c^3 + 37 c^4 + c^5) -
42 a^4 b^3 (-1 + c)^4 (1 + c)^3 -
30 a^4 b^2 (1 + c)^2 (11 + 76 c + 66 c^2 + 76 c^3 + 11 c^4) +
21 a b^6 (-1 + c)^4 (1 + c)^3 -
60 a b^5 (1 + c)^2 (7 - 4 c + 42 c^2 - 4 c^3 + 7 c^4) +
21 a^2 b^5 (-1 + c)^4 (1 + c)^3 -
30 a^2 b^4 (1 + c)^2 (11 + 76 c + 66 c^2 + 76 c^3 + 11 c^4) -
30 b^6 (1 + c)^2 (5 - 12 c + 30 c^2 - 12 c^3 + 5 c^4)) x -
a b (1 + c) (54 a^3 b^2 (-1 + c)^4 (1 + c) -
168 a^3 b (1 + c) (1 + 11 c + 11 c^2 + c^3) -
54 a^4 b (-1 + c)^4 (1 + c) -
6 a^4 (29 - 4 c + 286 c^2 - 4 c^3 + 29 c^4) -
54 a b^4 (-1 + c)^4 (1 + c) -
168 a b^3 (1 + c)^2 (1 + 10 c + c^2) +
54 a^2 b^3 (-1 + c)^4 (1 + c) +
12 a^2 b^2 (1 - 340 c - 330 c^2 - 340 c^3 + c^4) -
6 b^4 (29 - 4 c + 286 c^2 - 4 c^3 + 29 c^4)) x^2 -
a b (1 + c) (-20 a b (9 + 50 c + 9 c^2) -
10 a^2 (9 + 50 c + 9 c^2) - 10 b^2 (9 + 50 c + 9 c^2)) x^3 +
72 a b (1 + c) x^4

-
 Can you give us the polynomial explicitly..? – PlatoManiac Sep 16 '12 at 23:50 I did not at first because it is really ugly. – Apatura Sep 16 '12 at 23:55 "Number of roots" as in "number of real roots"? A cubic always has three roots, for instance... – J. M.♦ Sep 17 '12 at 0:02 Number of real roots >0. – Apatura Sep 17 '12 at 0:04

p[x_, a_, b_, c_] := the formula


I'd do this :

Manipulate[ CountRoots[ p[x, a, b, c], {x, 0, Infinity}],
{a, -100, 100}, {b, -100, 100}, {c, -100, 100}]


If a plot is needed one can proceed this way :

ListContourPlot3D[ Table[ CountRoots[ p[x, a, b, c], {x, 0, Infinity}],
{a, -10, 10}, {b, -10, 10}, {c, -10, 10} ], Contours -> 3,
Mesh -> None, ContourStyle -> {Red, Yellow, Lighter @ Blue},
DataRange -> {{-10, 10}, {-10, 10}, {-10, 10}}]


Another way (the most expensive) is to make use of RegionPlot3D. Here are 3d regions where a given polynomial has respectively : at least 1, 2 and 3 roots for x > 0

GraphicsRow[
Table[ RegionPlot3D[ CountRoots[ p[x, a, b, c], {x, 0, Infinity}] >= k,
{a, -10, 10}, {b, -10, 10}, {c, -10, 10},
PlotStyle -> Directive[ Orange, Opacity[0.5], Specularity[White, 30]]],
{k,  3}]]


One can see that for for higher k we need a better resolution, nevertheless it appears very expensive to compute regions using higher PlotPoints and MaxRecursion options.

 RegionPlot3D[ CountRoots[ p[x, a, b, c], {x, 0, Infinity}] >= 2,
{a, -10, 10}, {b, -10, 10}, {c, -10, 10},
PlotStyle -> Directive[ Orange, Opacity[0.5], Specularity[White, 30]],
PlotPoints -> 40, MaxRecursion -> 2]


-
 Thanks a lot, that is great for me to get to know the function! But I really need a plot I can print on paper. – Apatura Sep 17 '12 at 0:29 @Apatura I updated my answer with a ListContourPlot3D. – Artes Sep 17 '12 at 0:45 Cute plot work! – belisarius Sep 17 '12 at 4:53 @Artes: RegionPlot does exactly what I want, thank you very much! Unfortunately, when I try a different function, it gives back a strange error. What is wrong with RegionPlot3D[ CountRoots[\[Sigma]q33[\[Omega], V1, V2, V3], {\[Omega], 0, Infinity}] == 1, {V1, 0, 10}, {V2, 0, 10}, {V3, 0, 10}, PlotStyle -> Directive[Gray, Opacity[0.5]], Mesh -> None] if [Sigma]q33 is a polynomial in omega, Exp[-V1],Exp[-V2] and Exp[-V3]? – Apatura Sep 17 '12 at 14:52 @Apatura You are welcome. Note that when I set PlotPoints -> 60, MaxRecursion -> 3 it took more than 5 minutes to evaluate. If you have an error most likely it is because of a wrong definition of your function \[Sigma]q33. Look carefully at the documentation of RegionPlot3D. – Artes Sep 17 '12 at 14:57
show 1 more comment

Here's one approach. Generate data (for a simple polynomial):

data = Table[CountRoots[a*x^2 + b*x + c, x],
{a, -5, 5, 0.5}, {b, -5, 5, 0.5}, {c, -5, 5, 0.5}];


Then we'll display a collection of Cuboids with Opacity based on this data.

Graphics3D[{Opacity[0.2], EdgeForm[],
MapIndexed[Which[
# == 0, {Red, Cuboid[#2, #2 + {1, 1, 1}]},
# == 1, {White, Cuboid[#2, #2 + {1, 1, 1}]},
# == 2, {Blue, Cuboid[#2, #2 + {1, 1, 1}]}] &,
data, {3}]}]


Now it appears to turn out, that your more complicated functions all have zero, one, or two roots. Changing the With statement accordingly, we get the following.

You might also look into CUDAVolumetricRender, if you have a good NVidia graphics card.

-
 Might also be a good idea to generate a separate plot for the different root counts... – J. M.♦ Sep 17 '12 at 0:15 @J.M. I'm not sure I follow - that image is based on root counts. I turns out that all polynomials in that family have on zero, two, or four roots - hence, the three colors. – Mark McClure Sep 17 '12 at 0:19 Hmm, I thought you were doing something like ContourPlot3D[(* root-counting function *), {a, ...}, {b, ...}, {c, ...}, Contours -> {0, 2, 4}]... – J. M.♦ Sep 17 '12 at 0:22 I'd have used Switch[] instead of Which[] myself... something like {Switch[#, 0, Red, 1, White, 2, Blue], Cuboid[#2, #2 + {1, 1, 1}]} &. – J. M.♦ Sep 17 '12 at 0:39 @MarkMcClure will you please explain a bit how to use CUDAVolumetricRender in this context? – PlatoManiac Sep 17 '12 at 8:38