1
$\begingroup$

I'm working with generative polynomials, and would like to know if there is a faster way to calculate roots. Using the technique below I've managed to cut the time down significantly. However at n=35, it still takes about two and a half minutes to calculate roots. Does anyone know of any ways I could speed this up? The time per polynomial increases rather quickly, and I would like to be able to calculate up to n=50, and I do not have access to a very fast machine.

h = Compile[{{z, _Real}, {m, _Integer}}, 
Coefficient[Series[1/((1 + l t + z t^2)^2), {t, 0, m}], t^m]];
l = 6.0;
Do[{
roots = z /. {ToRules[Quiet[NRoots[h[z, n] == 0., z]]]},
Print[roots]
    }, {n, 2, 35}]
$\endgroup$
  • $\begingroup$ Do you need all roots? $\endgroup$ – Marius Ladegård Meyer May 3 '17 at 21:46
  • $\begingroup$ Take a look at First positive root $\endgroup$ – Artes May 4 '17 at 13:01
  • $\begingroup$ Could define h2[z_, n_] := Coefficient[ Sum[(-1)^k*Expand[(1 + l t + z t^2)^2 - 1]^k, {k, Ceiling[n/2], n}], t^n]. That avoids the Series bottleneck. There are several other ways to do this, as shown in various responses. $\endgroup$ – Daniel Lichtblau May 4 '17 at 15:46
3
$\begingroup$

Why not use SeriesCoefficient?

coeff[l_,n_] = Assuming[n>0, SeriesCoefficient[1/((1+l t+z t^2)^2),{t,0,n}]];
coeff[l,n] //TeXForm

$\frac{\left(\sqrt{l^2-4 z}+l\right)^2 \left(n \sqrt{l^2-4 z}+2 \sqrt{l^2-4 z}-l\right) \left(-\frac{1}{2} \sqrt{l^2-4 z}-\frac{l}{2}\right)^n}{4 \left(l^2-4 z\right)^{3/2}}+\frac{2^{-n-2} \left(\sqrt{l^2-4 z}-l\right)^{n+2} \left(n \sqrt{l^2-4 z}+2 \sqrt{l^2-4 z}+l\right)}{\left(l^2-4 z\right)^{3/2}}$

Solving the roots of coeff[l,n] is pretty fast:

roots = Table[
    NSolve@coeff[6, n],
    {n,2,50}
];// AbsoluteTiming
roots[[;;5]]

{2.6551, Null}

{{{z -> 54.}}, {{z -> 24.}}, {{z -> 126.991}, {z -> 17.0091}}, {{z -> 45.8745}, {z -> 14.1255}}, {{z -> 229.135}, {z -> 28.2533}, {z -> 12.612}}}

$\endgroup$
2
$\begingroup$

It looks like you are creating a series and then picking off the m-most term, which is a polynomial in z, and finding where it is equal to zero, correct? Just compute the m-most term directly from Taylor series properties:

hh[m_] := (D[1/((1 + l t + z t^2)^2), {t, m}] /. t -> 0)/m!

Then do the Do[ ] loop.

Do[{roots = z /. {ToRules[NRoots[hh[n] == 0., z]]}, Print[roots]}, {n,2, 50}]

Timing[ ] gave me .03 seconds to do all 50. Same answer as yours to do 25, but yours took 2.13 seconds just to do 25 of them.

$\endgroup$
  • $\begingroup$ Why use Timing[]? $\endgroup$ – Michael E2 May 4 '17 at 1:39
  • $\begingroup$ This was ridiculously fast, however it came out with a bunch of incorrect results at higher orders. Thanks though. $\endgroup$ – Daniel May 4 '17 at 1:41
  • $\begingroup$ @Daniel, the code I wrote is correctly computing the term for the Taylor Series, just doing it directly rather than through calling the series function. I compared the results of the terms out to n=35 for the uncompiled version of your function and my direct computation, and they are identical. Too slow to go to 50. If l is in fact an integer l=6 then there is no time penalty for leaving it that way, and you are avoiding floating point errors which may well be the source of difference between our codes. Keep your numbers exact as long as possible in your computations. $\endgroup$ – MikeY May 4 '17 at 12:53
0
$\begingroup$

Answer to the answer of Carl Woll:

There are two issus with that solution. 1. With MMA version 8.0 NSolve does't find all roots. Applying Simplify solves that problem.

Compare

 Table[{n, {NSolve[coeff[6, n]] // Sort}}, {n, 2, 
    10}] // TableForm

with

 Table[{n, {NSolve[coeff[6, n] // Simplify] // Sort}}, {n, 2, 
   10}] // TableForm
  1. Using NRoots instead of NSolve is up to 50 times faster. Compare

    (roots = 
       Table[{n, {NSolve[coeff[6, n] // Simplify] // Sort}}, {n, 2, 
          50}]); // Timing
    
    (*    {3.735, Null}   *)
    

with

  (roots2 = 
Table[{n, {{ToRules[NRoots[Simplify[coeff[6, n]] == 0, z]]} // 
    Sort}}, {n, 2, 50}]); // Timing

(*  {0.078, Null}    *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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