2
$\begingroup$

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 θ: 

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

Series expression of a Root object

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

ToRadicals@Roots

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?

Improve this question
$\endgroup$
0

4 Answers 4

Reset to default
2
$\begingroup$

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 *)

takes about an hour.

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.

Addendum: Numerical Expansion

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"}]

enter image description here

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – JI-br
    Commented Feb 16, 2017 at 9:02
  • $\begingroup$ @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, $\endgroup$
    – bbgodfrey
    Commented Feb 16, 2017 at 14:58
  • $\begingroup$ 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. $\endgroup$
    – JI-br
    Commented Feb 16, 2017 at 15:18
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Feb 17, 2017 at 6:08
  • $\begingroup$ Thank you @bbgodfrey. Your comprehensive explanation was of great help. $\endgroup$
    – JI-br
    Commented Mar 8, 2017 at 8:52
1
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Not a bad approximation to the symbolically derived answer. $\endgroup$
    – bbgodfrey
    Commented Feb 15, 2017 at 14:19
  • $\begingroup$ @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... $\endgroup$
    – JI-br
    Commented Feb 16, 2017 at 9:09
  • $\begingroup$ 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 $\endgroup$
    – JI-br
    Commented Feb 16, 2017 at 9:23
  • $\begingroup$ Tx for your help, I really appreciate that. $\endgroup$
    – JI-br
    Commented Mar 8, 2017 at 8:50
1
$\begingroup$

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}
Improve this answer
$\endgroup$
1
$\begingroup$

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.

Improve this answer
$\endgroup$

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.