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
$\endgroup$
4
  • $\begingroup$ Can you give us the polynomial explicitly..? $\endgroup$ Commented Sep 16, 2012 at 23:50
  • $\begingroup$ I did not at first because it is really ugly. $\endgroup$
    – Apatura
    Commented 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$ Commented Sep 17, 2012 at 0:02
  • $\begingroup$ Number of real roots >0. $\endgroup$
    – Apatura
    Commented Sep 17, 2012 at 0:04

2 Answers 2

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

$\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
    Commented Sep 17, 2012 at 0:29
  • $\begingroup$ @Apatura I updated my answer with a ListContourPlot3D. $\endgroup$
    – Artes
    Commented Sep 17, 2012 at 0:45
  • $\begingroup$ Cute plot work! $\endgroup$ Commented 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
    Commented 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
    Commented 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.

$\endgroup$
7
  • $\begingroup$ Might also be a good idea to generate a separate plot for the different root counts... $\endgroup$ Commented 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$ Commented 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$ Commented 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$ Commented Sep 17, 2012 at 0:39
  • $\begingroup$ @MarkMcClure will you please explain a bit how to use CUDAVolumetricRender in this context? $\endgroup$ Commented Sep 17, 2012 at 8:38

Your Answer

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

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