Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

How to solve this integral by Mathematica even by numerical methods (plotting the solution)

Integrate[(Cos[x] - a)/(1 + a^2 - 2*a*Cos[x])^1.5, {x, 0, 2*Pi}]

enter image description here

It consumes huge time on Mathematica 9 home edition and no output.

I just want the Taylor expansion in the vicinity of a=0

share|improve this question
    
Why integral over [0,1] when you have [0, 2 Pi] in the picture ? –  b.gatessucks Oct 23 '13 at 20:24
    
Yes, it was an error. I edited the post to be [0, 2 Pi]. Thanks –  أحمد كمال Oct 24 '13 at 2:08

3 Answers 3

expr =  Assuming[a > 0, Integrate[(Cos[x] - a)/(1 + a^2 - 2*a*Cos[x])^(3/2), {x, 0, 1}]]

gives:

(-(-1 + a) Sqrt[1 + a^2 - 2 a Cos[1]] ((-1 + a) EllipticE[1/
      2, -((4 a)/(-1 + a)^2)] + (1 + a) EllipticF[1/
      2, -((4 a)/(-1 + a)^2)]) - 
 2 a Abs[-1 + a] Sin[1])/(a (-1 + a^2) Abs[-1 + a] Sqrt[
 1 + a^2 - 2 a Cos[1]])

Note: I replaced 1.5 by the radical 3/2. If you want numerical value use NIntegrate

As suggested by b.gatessuck (also belisarius' answer), taking the indefinite integral makes this faster:

expr = Integrate[(Cos[x] - a)/(1 + a^2 - 2*a*Cos[x])^(3/2), x]

Gives:

 (-(-1 + a)^2 Sqrt[(1 + a^2 - 2 a Cos[x])/(-1 + a)^2]
       EllipticE[x/2, -((4 a)/(-1 + a)^2)] - (-1 + a^2) Sqrt[(
      1 + a^2 - 2 a Cos[x])/(-1 + a)^2]
       EllipticF[x/2, -((4 a)/(-1 + a)^2)] - 
     2 a Sin[x])/(a (-1 + a^2) Sqrt[1 + a^2 - 2 a Cos[x]])

You can then Plot this like any other function:

Plot[(expr/.x -> 1) - ( expr /. x -> 0 ), {a, 0, 5}]

Mathematica graphics

OR from 0 to 2 Pi

Plot[(expr/.x -> 2 Pi) - ( expr /. x -> 0 ), {a, 0, 5}]

Mathematica graphics

share|improve this answer
    
Try the indefinite integral first, then take the limits and see if it's faster. –  b.gatessucks Oct 23 '13 at 20:25
    
@b.gatessucks. Thanks, that definitely makes it instantaneous. –  RunnyKine Oct 23 '13 at 20:30
    
And that was my answer :) –  belisarius Oct 23 '13 at 20:32
    
@belisarius Beat you by a whole 35 seconds. –  b.gatessucks Oct 23 '13 at 20:32
1  
@b.gatessucks We should take some relativistic effects in account. Here G is fluctuating madly right now –  belisarius Oct 23 '13 at 20:34

This is instantaneous:

p = Integrate[(Cos[x] - a)/(1 + a^2 - 2*a*Cos[x])^(3/2), x, 
    Assumptions -> Element[a, Reals], GenerateConditions -> True];
Plot[(p /. x -> 1) - (p /. x -> 0), {a, 0, 5}]

Mathematica graphics

In case you made a mistake and the upper integration limit was 2 Pi:

Plot[(p /. x -> 2 Pi) - (p /. x -> 0), {a, 0, 5}]

Mathematica graphics

Edit

As for the series expansion:

s = Normal@Series[(p /. x -> 2 Pi) - (p /. x -> 0), {a, 0, 10}]

a π + (9 a^3 π)/8 + (75 a^5 π)/64 + (1225 a^7 π)/1024 + (19845 a^9 π)/16384

The numerators are the Denominators of the column sums of the BG2 matrix

The denominators follows some power of 2 rule, but not trivial:

s = Series[(p /. x -> 2 Pi) - (p /. x -> 0), {a, 0, 100}];
lf = Denominator /@ List @@ Normal@s;
ListLogPlot@lf

Mathematica graphics

Perhaps the following may help you if you need to find a closed expression:

s = Series[(p /. x -> 2 Pi) - (p /. x -> 0), {a, 0, 1000}];
lf = Differences[Last /@ Flatten /@ FactorInteger /@ Denominator /@ List @@ Normal@s];
ListLinePlot@lf

Mathematica graphics

share|improve this answer
1  
Just p = Integrate[(Cos[x] - a)/(1 + a^2 - 2*a*Cos[x])^(3/2), x] will do, other options seem to be redundant in this case. –  Rorschach Oct 23 '13 at 20:37
1  
@Blackbird GenerateConditions is a safeguard against special values, so I'd rather prefer to let it be. The Reals assumption is just a clichè –  belisarius Oct 23 '13 at 20:41

This will compute the Taylor expansion:

Define the integrand.

g[θ_, a_] := (Cos[θ] - a)/(1 + a^2 - 2 a Cos[θ])^(3/2)

Integrate a series expansion.

gi[a_, n_] := Integrate[Series[g[θ, a], {a, 0, n}] // Normal, {θ, 0, 2 π}]

Find a sequence function that fits a sufficiently long series expansion.

gs[a_, n_] = FindSequenceFunction[List @@ Expand@gi[a, 15]][n] // FullSimplify

(* (4 (a^2)^n Gamma[1/2 + n]^2)/(a Gamma[n] Gamma[1 + n]) *)

Summing this over $n = 1, 2, ...$ gives you the required series expansion.

share|improve this answer
    
Nice one :) I tried FindSequenceFunction[] with 12 terms and then abandoned the idea. –  belisarius Oct 24 '13 at 3:09

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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