5
$\begingroup$

Here is my problem ; for my research, I believe that the complex numbers I am looking at are precisely the (very large) set of roots of some high degree polynomial, of degree $\sim 2^n$ where $1 \le n \le 10 \sim 15$. Mathematica has been running for the whole day on my computer just for $2^{10}$ even though $2^9$ took half an hour, and I wondered if any other program out there would be faster than Mathematica so that I could compute more of those roots. If I had more examples to compute it would REALLY help. The thing I need the program to be able to do is simple : I give you a polynomial of very high degree, and I want to numerically plot its complex roots. I don't care about multiplicity.

Thanks in advance,

EDIT: Marvis, here is my code for computing the nested polynomial $p^m(x) - x$, where $p^2(x) = p(p(x))$.

function r = polycomp(p,q);

r = p(1);

for k = 2:length(p);

r = conv(r,q);

r(end) = r(end) + p(k);

end

All I do afterwards is a loop with

r = [1 0]

for i = 1:n

r = polycomp(r,p)

end

where $n$ is my loop length and $p$ is my polynomial.

  • 0
    Have you tried using MATLAB? If not, try and use the command roots http://www.mathworks.com/help/techdoc/ref/roots.html. I tried for polynomial of degree 2000 and it took around a minute on my computer.2012-06-23
  • 0
    @Marvis : Seriously? Man that is WAY faster than Mathematica! I didn't expect MATLAB would be so good.2012-06-23
  • 0
    @Marvis : Do you know any command in MATLAB allowing the "nesting" of polynomials? I don't expect this is a computational problem but to compute my polynomial I know it has the form $p(p(p(p(...(x))...)$ for a quadratic $p$, so it is unlikely I actually know the coefficients myself of the polynomial I am looking for..2012-06-23
  • 1
    I have not solved polynomials in MATLAB so I am not much aware of solving polynomials. But I was wondering if you wanted to solve say $p(p(x)) = 0$, won't it be better to solve $p(y^*) = 0$ and then solve $p(x) = y^*$? You can then do this recursively. If you do it this way, you do not need to compute the coefficients. This is just a suggestion and I am not completely sure if this a better way to go.2012-06-23
  • 0
    Actually I am solving $p(p(...p(x)...)) = x$, so no, a nesting is really not to be expected.2012-06-23
  • 0
    @PatrickDaSilva MATLAB represents polynomials as vectors of coefficients, so it should be easy to write a simple *for* loop that computes the coefficients iteratively. Good MATLAB exercise, too.2012-06-23
  • 0
    @Marvis : I see that MATLAB is really fast, but I have problems with precision ; at some point the coefficients of my polynomial become "INF" because they are too big... is there anyway I could ask matlab to raise its precision?2012-06-23
  • 0
    @Leonid : I figured it out, I just wondered if there was some neat function. I made one myself using *conv*, it's okay.2012-06-23
  • 0
    @PatrickDaSilva Oh yes this will be an issue unfortunately. Is there a way to keep the maximum of the absolute value of your coefficient to be $1$? At the moment, I can only think of this.2012-06-23
  • 0
    Well I only want to plot examples, so I could choose "computable" examples, but won't the opposite problem happen if I choose coefficients too small, i.e. they'll run to zero?2012-06-23
  • 0
    @Marvis : Actually I managed to do it for $2^{11}$! This is a first. Just for this, thanks for the tips! I chose the coefficients of my quadratic smaller than $1$ in absolute value and it worked. It gives me a beautiful plot! If you want to post the main point of your arguments in an answer I'd be willing to +1/check it2012-06-23
  • 0
    I don't remember the code to put everything in a gray box as above... it magically did itself and I can't put everything in it. Anyone knows the command? What I have right now is plenty of one line boxes. >.>2012-06-23
  • 0
    Perhaps you could also try something like this, where you plot the whole function in such a way that it's visually clear where the zeros are: http://www.mai.liu.se/~halun/complex/domain_coloring-unicode.html#fig:iterates2012-06-24

2 Answers 2

2

The MATLAB command roots() seems to do a fine job and at a decent speed. I tried using roots for polynomials of degree from $2^7$ to $2^{14}$ and below are the timings. The time for computing the roots using roots() seem to scale as $N^{\log_2(5)}$. enter image description hereenter image description here

  • 0
    Thanks for everything. But now I am facing a new problem to go above $2^{13}$ (although I am VERY HAPPY to have gotten there! I have beautiful graphs I can conjecture at now.) The problem is I can't compute, say Nest[p,x,14] - x precisely enough. If the coefficients are integers, they go to MATLAB's infinity : if they are all $|a|,|b|,|c|< 1$ for say $ax^2 + bx + c$, then they go to zero too fast, and become NaN's. I have approximately your specs, but I am not worried about the time it takes, because I can wait (as long as its reasonable time, say a day or less... I only have one life left!)2012-06-23
  • 0
    How do you compute the coefficients of the final polynomial? I am wondering if there is any trick we can do to avoid computing these large coefficients in the first place...2012-06-23
  • 0
    Well I can show you my code if you want! Comments are bad for this. See my answer's edit2012-06-23
1

I'm coming to the conversation late, but I note that there is no discussion of the accuracy of the determined roots. Roots are highly sensitive to perturbations in the polynomial coefficients, i.e. numerical precision is a key determinant of root accuracy particularly for huge polynomials.

My question: MATLAB may be faster but is it more accurate than Mathematica? Did you check to see whether your $2^n$ roots aren't mostly wildly inaccurate?

  • 0
    I am not working on the problem relative to this question anymore, and I managed to compute the roots to sufficient precision to notice that the problem I was working on is not of interest. Still, thank you for answering.2013-01-19