6
$\begingroup$

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
Improve this question
$\endgroup$
4
  • $\begingroup$ Can you give us the polynomial explicitly..? $\endgroup$
    – PlatoManiac
    Sep 16, 2012 at 23:50
  • $\begingroup$ I did not at first because it is really ugly. $\endgroup$
    – Apatura
    Sep 16, 2012 at 23:55
  • $\begingroup$ "Number of roots" as in "number of real roots"? A cubic always has three roots, for instance... $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 17, 2012 at 0:02
  • $\begingroup$ Number of real roots >0. $\endgroup$
    – Apatura
    Sep 17, 2012 at 0:04

2 Answers 2

Reset to default
10
$\begingroup$

Defining your polynomial as 

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}]

enter image description here

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}}]

enter image description here

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}]]

enter image description here

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]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Apatura
    Sep 17, 2012 at 0:29
  • $\begingroup$ @Apatura I updated my answer with a ListContourPlot3D. $\endgroup$
    – Artes
    Sep 17, 2012 at 0:45
  • $\begingroup$ Cute plot work! $\endgroup$
    – Dr. belisarius
    Sep 17, 2012 at 4:53
  • $\begingroup$ @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]? $\endgroup$
    – Apatura
    Sep 17, 2012 at 14:52
  • $\begingroup$ @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. $\endgroup$
    – Artes
    Sep 17, 2012 at 14:57
7
$\begingroup$

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}]}]

enter image description here

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.

enter image description here

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

Improve this answer
$\endgroup$
7
  • $\begingroup$ Might also be a good idea to generate a separate plot for the different root counts... $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 17, 2012 at 0:15
  • $\begingroup$ @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. $\endgroup$
    – Mark McClure
    Sep 17, 2012 at 0:19
  • $\begingroup$ Hmm, I thought you were doing something like ContourPlot3D[(* root-counting function *), {a, ...}, {b, ...}, {c, ...}, Contours -> {0, 2, 4}]... $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 17, 2012 at 0:22
  • $\begingroup$ I'd have used Switch[] instead of Which[] myself... something like {Switch[#, 0, Red, 1, White, 2, Blue], Cuboid[#2, #2 + {1, 1, 1}]} &. $\endgroup$
    – J. M.'s eventual burnout ♦
    Sep 17, 2012 at 0:39
  • $\begingroup$ @MarkMcClure will you please explain a bit how to use CUDAVolumetricRender in this context? $\endgroup$
    – PlatoManiac
    Sep 17, 2012 at 8:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.