# Series expression of a RootSum object

Mathematica was not able to calculate the definite integral of a trigonometric function of mine:

Assuming[r0 > 0,
Assuming[r0 < 1/2,
Integrate[(Cos[θ]^6 Sin[θ]^2)/(a + Cos[θ]^3), {θ, 0, 2 π}]]]


so I made a standard change of variables Tan[θ/2]=z, transforming the trigonometric function in a rational function in the variable z. Then I calculate the indefinite integral in z (I also replace r0 by a, using the relation a -> (-1 + 1/r0^3) )

I1 = Integrate[(4 z^2 (1 - z^2)^6)/((1 + z^2)^8 (a + (1 - z^2)^3/(1 + z^2)^3))*2 /(1 + z^2), z]


Mathematica returned me the calculation in terms of RootSum

1/15 ((-15 a z (-1 + z^2) (1 + z^2)^3 + 8 z^3 (5 - 2 z^2 + 5 z^4) -
15 a (1 + z^2)^5 ArcTan[z])/(1 + z^2)^5 +
20 a^2 RootSum[
1 + a - 3 #1^2 + 3 a #1^2 + 3 #1^4 + 3 a #1^4 - #1^6 +
a #1^6 &, (
Log[z - #1] #1)/(-1 + a + 2 #1^2 + 2 a #1^2 - #1^4 + a #1^4) &])


and then I made the way back so I would have the result in terms of the original variable Theta:

I2 = I1 /. ArcTan[z] -> θ/2;

I3 = I2 /. z -> Sin[θ]/(1 + Cos[θ])

1/15 (20 a^2 RootSum[
1 + a - 3 #1^2 + 3 a #1^2 + 3 #1^4 + 3 a #1^4 - #1^6 +
a #1^6 &, (
Log[Sin[θ]/(1 + Cos[θ]) - #1] #1)/(-1 + a +
2 #1^2 + 2 a #1^2 - #1^4 + a #1^4) &] + (-((
15 a Sin[θ] (-1 +
Sin[θ]^2/(1 + Cos[θ])^2) (1 +
Sin[θ]^2/(1 + Cos[θ])^2)^3)/(
1 + Cos[θ])) -
15/2 a θ (1 + Sin[θ]^2/(1 + Cos[θ])^2)^5 + (
8 Sin[θ]^3 (5 - (
2 Sin[θ]^2)/(1 + Cos[θ])^2 + (
5 Sin[θ]^4)/(1 + Cos[θ])^4))/(1 +
Cos[θ])^3)/(1 + Sin[θ]^2/(1 + Cos[θ])^2)^5)


which I simplified using FullSimplify

I33 = FullSimplify[I3]

1/240 (40 a (-3 θ +
8 a RootSum[
1 + a - 3 #1^2 + 3 a #1^2 + 3 #1^4 + 3 a #1^4 - #1^6 +
a #1^6 &, (
Log[-#1 + Tan[θ/2]] #1)/(-1 + a + 2 #1^2 +
2 a #1^2 - #1^4 + a #1^4) &]) + 30 Sin[θ] +
60 a Sin[2 θ] - 5 Sin[3 θ] - 3 Sin[5 θ])


Finally I calculated the integral between 0 and 2π

II = FullSimplify[
2 (Limit[I33, θ -> π,
Direction -> 1] - (I33 /. θ -> 0))]

-(1/3)
a (3 π +
8 a RootSum[
1 + a - 3 #1^2 + 3 a #1^2 + 3 #1^4 + 3 a #1^4 - #1^6 +
a #1^6 &, (
Log[-#1] #1)/(-1 + a + 2 #1^2 + 2 a #1^2 - #1^4 + a #1^4) &]);

IIFINAL /. a -> (-1 + 1/r0^3)

-(1/
3) (-1 + 1/r0^3) (3 π -
8 (-1 + 1/
r0^3) r0^3 RootSum[-1 - 3 #1^2 + 6 r0^3 #1^2 - 3 #1^4 - #1^6 +
2 r0^3 #1^6 &, (
Log[-#1] #1)/(-1 + 2 r0^3 - 2 #1^2 - #1^4 + 2 r0^3 #1^4) &])


Now I want to calculate the Series expansion of this expression, in terms of r0.

I would like to do something like

Series[IIFINAL, {r0, 0, 5}]


but when I try to do that, Mathematica returns me

Series::nmer: "Root[-1+(-3+6\ r0^3)\ #1-3\ #1^2+(-1+2\ r0^3)\ #1^3&,1] is not a meromorphic function of r0 at 0."

I have already checked the post

but it does not help me since in my case I have RootSum, not Root in the result, and when I try to use

as suggested in that post, I get an error

Nonpolynomial root specification .... is not a list with two elements.

So, is there a way to force Mathematica to output a series representation for RootSum-like results?

Symbolic Expansion

With the final expression in the question designated as f, a series expansion can be obtained by converting the RootSum to a set of Roots, converting them ToRadicals, generating a Series, and converting its coefficients to numerical values. Results of the intermediate steps are too lengthy to reproduce here. The entire code is

Series[ToRadicals@Normal@f, {r0, 0, 10}] // Normal // N // Chop


which after several minutes returns

(* 0.245437 r0^3 + 0.245437 r0^6 + 0.34668 r0^9 *)


Obtaining another term,

Series[ToRadicals@Normal@f, {r0, 0, 13}] // Normal // N // Chop
(* 0.245437 r0^3 + 0.245437 r0^6 + 0.34668 r0^9 + 0.549165 r0^12 *)


Note that these Series contain only powers of r0^3, consistent with the fact that r0 appears in f only in terms of r0^3. This suggests replacing r0^3 by z.

Series[ToRadicals@Normal@(f /. r0 -> z^(1/3)), {z, 0, 4}] // Normal // N // Chop
(* 0.245437 z + 0.245437 z^2 + 0.34668 z^3 + 0.549165 z^4 *)


Unfortunately, this transformation yields no improvement in speed.

A numerical expansion yields results much more rapidly.

<< NumericalCalculus
Chop[NSeries[f /. r0 -> z^(1/3), {z, 0, 20}, Radius -> .3,
WorkingPrecision -> 60], 10^-4] // Normal
(* 0.245437 z + 0.245437 z^2 + 0.34668 z^3 + 0.549165 z^4 + 0.911161 z^5 +
1.5492 z^6 + 2.67704 z^7 + 4.68355 z^8 + 8.27794 z^9 + 14.759 z^10 +
26.5155 z^11 + 47.9579 z^12 + 87.2601 z^13 + 159.622 z^14 + 293.398 z^15 +
541.633 z^16 + 1003.83 z^17 + 1867.12 z^18 + 3484.19 z^19 + 6521.26 z^20 *)


Because the Series coefficients increase steadily, it seems prudent to test convergence.

Plot[{Re@Evaluate[fun /. r0 -> z^(1/3)], %}, {z, -1/2, 1/2}, PlotRange -> All,
AxesLabel -> {"z", "f"}]


Agreement is very good except near z == 0.5, where a branch point occurs.

• Thank you @bbgodfrey. I tried to follow your steps but Mathematica is not actually able to calculate the Series part. I mean, I kept it running all night long and it was not able to deliver an answer. My computer is a notebook Dell with a 2.66GHz Intel i5 processor with 8GB of RAM running Wofram Mathematica 9. Should I need more RAM to run this? – JI-br Feb 16 '17 at 9:02
• @JI-br The calculation is somewhat slow but does not need a lot of RAM. However, I am using Mathematica 11.0.1 on Windows 10 (64-bit), and you are using one of the version 9's, which may well make a difference. By the way, the series must be in powers of r0^3, and I am trying a slightly different approach right now that takes advantage of that fact, – bbgodfrey Feb 16 '17 at 14:58
• I understand that previous versions of Mathematica might yeld additional complication in the calculations. Now, how do you know that the series must be in powers of r0^3? Using the strategy proposed by @Manu but expanding the series up to the 20th power, it seems that all the powers of r0 greater or equal 9 are significant in the series. – JI-br Feb 16 '17 at 15:18
• @JI-br In your final expression (the one you wish to expand in a series), r0 appears only as powers of r0^3, so the same must be true for the series. The approximation employed by Manu breaks down at larger terms, although increased WorkingPrecision might help. – bbgodfrey Feb 17 '17 at 6:08
• Thank you @bbgodfrey. Your comprehensive explanation was of great help. – JI-br Mar 8 '17 at 8:52

I am not an expert, so someone else is waiting for an authoritative answer.

I just wanted to make you see that writing:

f[r_?NumericQ] := NIntegrate[Cos[t]^6 Sin[t]^2 / (Cos[t]^3 - 1 + 1/r^3), {t, 0, 2 Pi}];

Chop[Series[f[r], {r, $MachineEpsilon, 9}], 10^-1]  I get: 0.245437 r^3 + 0.245463 r^6 + 0.346255 r^9 + O[r]^10  See if it is what you want. • Not a bad approximation to the symbolically derived answer. – bbgodfrey Feb 15 '17 at 14:19 • @Manu thank you for your answer. This solution runs fast in my computer. I wonder how much I lose by replacing zero to$MachineEpsilon, though... – JI-br Feb 16 '17 at 9:09
• One more question Manu and @bbgodfrey: I don't see why the computer quickly computes the answer when I do f[r_?NumericQ]= := NIntegrate[Cos[t]^6 Sin[t]^2 / (Cos[t]^3 - 1 + 1/r^3), {t, 0, 2 Pi}] Series[f[r], {r, $MachineEpsilon, 10}] // Chop but it is not able to provide an anwser (it keeps running forever) when I do f[r_]= := Integrate[Cos[t]^6 Sin[t]^2 / (Cos[t]^3 - 1 + 1/r^3), {t, 0, 2 Pi}] Series[f[r], {r,$MachineEpsilon, 10}] // Chop – JI-br Feb 16 '17 at 9:23
• Tx for your help, I really appreciate that. – JI-br Mar 8 '17 at 8:50

Here is a different way to go about expanding your integral into a series. Through symmetry considerations, your integral is equivalent to

2 Integrate[(Cos[θ]^6 Sin[θ]^2)/(1/r0^3 - 1 - Cos[θ]^3), {θ, 0, π}]


Now, we shall swap summation and integration (whose justification I will leave for someone else), and try to expand the integrand as a series in r0. Because of how it is constructed, we expect that the only nonzero terms would be the coefficients of $$x^{3k}$$:

Assuming[k > 0 && k ∈ Integers && -1 < c < 1 && -1 < s < 1,
FullSimplify[Simplify[
SeriesCoefficient[2 (c^6 s^2)/(1/r0^3 - 1 - c^3), {r0, 0, k}]] /.
k -> 3 k]] /. Thread[{c, s} -> {Cos[θ], Sin[θ]}]
2 Cos[θ]^6 (1 + Cos[θ]^3)^(-1 + k) Sin[θ]^2


Thus, one can quickly derive the following series:

With[{n = 6},
Sum[NIntegrate[2 Cos[θ]^6 (1 + Cos[θ]^3)^(k - 1) Sin[θ]^2, {θ, 0, π},
WorkingPrecision -> 20] x^(3 k), {k, 1, n},
Method -> "Procedural"] + O[x]^(3 n + 1)]
0.24543692606170259676 x^3 + 0.24543692606170259676 x^6 +
0.34667965806215491792 x^9 + 0.54916512206305956024 x^12 +
0.91116061955426017736 x^15 + 1.5492007535154440762 x^18 + O[x]^19


If one wants a symbolic answer, we need to do a little bit more. Recall that Sin[θ]^2 == 1 - Cos[θ]^2, and split the integral into two simpler ones:

Assuming[k >= 1,
2 Integrate[Cos[θ]^6 (1 + Cos[θ]^3)^(k - 1), {θ, 0, π}] -
2 Integrate[Cos[θ]^8 (1 + Cos[θ]^3)^(k - 1), {θ, 0, π}]]
5/8 π HypergeometricPFQ[{7/6, 3/2, 11/6, 1/2 - k/2, 1 - k/2},
{1/2, 4/3, 5/3, 2}, 1] -
35/64 π HypergeometricPFQ[{3/2, 11/6, 13/6, 1/2 - k/2, 1 - k/2},
{1/2, 5/3, 2, 7/3}, 1]


where we obtain an answer involving $${}_5 F_4$$ hypergeometric functions.

The hypergeometric expression itself does not seem to have a simpler form, but it surprisingly evaluates to more familiar expressions:

Table[5/8 π HypergeometricPFQ[{7/6, 3/2, 11/6, 1/2 - k/2, 1 - k/2},
{1/2, 4/3, 5/3, 2}, 1] -
35/64 π HypergeometricPFQ[{3/2, 11/6, 13/6, 1/2 - k/2, 1 - k/2},
{1/2, 5/3, 2, 7/3}, 1], {k, 9}]
{5 π/64, 5 π/64, 113 π/1024, 179 π/1024, 38015 π/131072, 64635 π/131072,
3574083 π/4194304, 6252949 π/4194304, 2829255293 π/1073741824}

N[%, 20]
{0.24543692606170259675, 0.24543692606170259675, 0.34667965806215491792,
0.54916512206305956024, 0.91116061955426017731, 1.5492007535154440763,
2.6770384064007208317, 4.6835466961077795126, 8.2779374379840101774}


While for speed reasons, JM's answer is certainly the way to go for higher order terms, if you are satisfied with the leading order terms you can use AsymptoticIntegrate:

AsymptoticIntegrate[
(Cos[θ]^6 Sin[θ]^2)/(1/r0^3-1-Cos[θ]^3),
{θ, 0, 2π},
{r0, 0, 10}
]
`

(5 π r0^3)/64 + (5 π r0^6)/64 + (113 π r0^9)/1024

reproducing JM's first 3 terms.