I have a simple code to get a planet's true anomaly as a function of time, using Newton's method and a series expansion. (Note that I purposefully iterate Newton's method two more times than the order I end up taking the expansion in)

I have used this code for years on M9 and 10. I just got 11, and now takes an obscenely long time to expand anything past linear.

    TA[k_]:= Block[{Qfunc, EoM, foE, foM},
        Qfunc[M_, e_, Ecc_] := -(Ecc - e Sin[Ecc] - M)/(1 - e Cos[Ecc]);

        EoM[e_, M_, 1] := M;
        EoM[e_, M_, kk_] := EoM[e, M, kk] = EoM[e, M, kk - 1] + Qfunc[M, e, EoM[e, M, kk - 1]];

        foE[e_, Ecc_] := 2 ArcTan[Sqrt[(1 + e)/(1 - e)] Tan[Ecc/2]] + Ecc - 2 ArcTan[Tan[Ecc/2]];
        foM[e_, M_, kk_] := foE[e, EoM[e, M, kk]];

        Series[foM[e, t n, k + 2], {e, 0, k}]


Then you run it with


The first two output quite quickly, the last one takes a very long time. On version 11, it looks like

{0.195877, Null}

{0.15612, Null}

{174.727, Null}

Why did it take so much longer to find the e^2 term? On previous versions of Mathematica with this exact same code it has always taken fractions of a second to expand even to much higher order.

Maybe there's an oversight, or some setting I need to tweak for this new version?

Thanks in advance!

  • $\begingroup$ Computing foM[e, t n, 4] takes negligible time. However, expanding this quantity to second order in e is very slow. $\endgroup$ – bbgodfrey Oct 11 '18 at 23:26
  • 1
    $\begingroup$ Seems to be an issue in version 11.3. It's fast again for the next release though. $\endgroup$ – Daniel Lichtblau Oct 22 '18 at 23:16
  • $\begingroup$ I ended up just updating to M12. It works now. $\endgroup$ – M.Walker May 9 '19 at 15:30

A bit of experimentation reveals that essentially all the time goes to computing Series rather than its arguments. Probably, Series spends this inordinate amount of time checking for branch cut issues. Using Assumptions -> (-1/2 < e < 1/2 && t n ∈ Reals) does not help. The following work-around produces fast results. Define

series[f_, {e_, e0_, n_}] := Sum[Simplify[D[f, {e, i}] /. e -> e0] e^i/i!, {i, 0, n}]

and use it instead of Series in the code in the question. Runtime is reduced by more than three orders of magnitude.

TA[2] // Timing
(* {0.078125, n t + 2 e Sin[n t] + 
              5/2 e^2 Cos[(n t)/2] (-Sin[(n t)/2] + Sin[(3 n t)/2])} *)

(Be sure to run Simplify on the answer obtained with Series before comparing results.)

| improve this answer | |

Rather than using Series on the whole expression, you can replace e with a SeriesData expression:

res = foM[e + O[e]^3, t n, 4] //Simplify; //AbsoluteTiming
res //TeXForm

{0.000169, Null}

$n t+2 e \sin (n t)+\frac{5}{2} e^2 \left(\sin \left(\frac{3 n t}{2}\right)-\sin \left(\frac{n t}{2}\right)\right) \cos \left(\frac{n t}{2}\right)+O\left(e^3\right)$

| improve this answer | |

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.