Re: (Not trivial) Definite Integration of a rational function
- To: mathgroup at smc.vnet.net
- Subject: [mg74811] Re: (Not trivial) Definite Integration of a rational function
- From: "dimitris" <dimmechan at yahoo.com>
- Date: Fri, 6 Apr 2007 04:20:12 -0400 (EDT)
- References: <euslob$9nf$1@smc.vnet.net>
Thanks a lot for anyone's rersponse. I think that when Mathematica integration algorithm deals with an improper integral does not use the Newton-Leibniz formula. Instead as the implementation notes state the integration proceeds using Marichev-Adamchik Mellin transform methods. The results are often initially expressed in terms of Meijer G functions, which are converted into hypergeometric functions using Slater's Theorem and then simplified. ( See http://library.wolfram.com/infocenter/Conferences/4684/ and also http://library.wolfram.com/infocenter/Conferences/5832/ ). So having in mind this fact I believe it can be explained why even though the indefinite integral is obtained the definite integral from zero to infinity returns unevalueted. I hope someone confirms or contradicts my belief! Dimitris =CF/=C7 dimitris =DD=E3=F1=E1=F8=E5: > I think it is time to open a new thread about integration, isn't it? > > Anyway...Here we go..(Mathematica 5.2 is used). > > > f = HoldForm[Integrate[(2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1)), {x, 0, > Infinity}]] > > We have > > ReleaseHold[f /. Integrate[h_, {x_, a_, b_}] :> Plot[h, {x, a, 10}]] > (*plot*) > > ReleaseHold[f /. Integrate -> NIntegrate] (*numerical estimation with > machine precision*) > 0.3712169752602472 > > Attempts to get an analytic result failed. > > Timing[ReleaseHold[f]] > {145.969*Second, Integrate[(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)), > {x, 0, Infinity}]} > > Timing[ReleaseHold[f /. Integrate[h___] :> Integrate[h, > GenerateConditions -> False]]] > {144.765*Second, Integrate[(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)), > {x, 0, Infinity}, GenerateConditions -> False]} > > ( > > BTW, Mathematica version 4 returns Infinity for this definite > integral. I think this happens because a partial fraction > decomposition is performed and the integrand is written > > Apart[(2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1))] > -(2/(1 + x)) + (2*(1 + 2*x + x^2))/(1 + 2*x + 3*x^2 + x^3) > > and after it integrates definitely each term of last sum > > (Integrate[#1, {x, 0, Infinity}] & ) /@ % > Integrate::idiv: Integral of 1/(1 + x) does not converge on > {0,Infinity}. > Integrate::idiv: Integral of (1 + x)^2/(1 + 2*x + 3*x^2 + x^3) does > not \ > converge on {0,Infinity}. > Integrate[-(2/(1 + x)), {x, 0, Infinity}] + Integrate[(2*(1 + 2*x + > x^2))/(1 + 2*x + 3*x^2 + x^3), {x, 0, Infinity}] > > One additional curious thing is it doesn't show the relevant warning > message for convergence... > ) > > Let see the antiderivative > > ReleaseHold[f /. Integrate[h_, {x_, a_, b_}] :> Integrate[h, x]] > Together[D[%, x]] > Plot[%%, {x, 0, 10}]; (*plot*) > > 2*(-Log[1 + x] + RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[x - #1] + > 2*Log[x - #1]*#1 + Log[x - #1]*#1^2)/ > (2 + 6*#1 + 3*#1^2) & ]) > (2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)) > > I wonder why Mathematica since it gets the indefinite integral, it > fails to evaluate > the definite integral. Any ideas? > > Since both the integrand and the antiderivative posses no > singularities/discontinuities > in the integration range I proceed to get the definite integral by > application of the Newton-Leibniz > formula > > This time we get finally the definite integral > > Timing[Limit[F, x -> Infinity] - Limit[F, x -> 0]] > Simplify[%[[2]]] > N[%]//Chop > > {176.90699999999998*Second, -2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , > (Log[-#1] + 2*Log[-#1]*#1 + Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) > & ]} > > -2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] + 2*Log[-#1]*#1 + > Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ] > 0.37121697526024766 > > But timings like this look quite unreasonable for integrals like this > (at least for me!). > > Do you have other ideas of getting this definite integral (with better > time performances > if possible)? > > It is also looks quite strange that in two other CAS (the one can be > freely downloaded from Internet) > I got directly the definite integral in a couple of seconds. > > Any explanation about previous time performance of Limit will be > greatly appreciate. > > Trying to write down the result > > -2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] + 2*Log[-#1]*#1 + > Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ] > > in a expanded format (i e no RootSum object; of course I am aware that > it will be lost the compacteness > that RootSum offers) I tried > > Tr[Simplify[(-2*((Log[-#1] + 2*Log[-#1]*#1 + Log[-#1]*#1^2)/(2 + 6*#1 > + 3*#1^2)) & ) /@ > (x /. Solve[(1 + 2*#1 + 3*#1^2 + #1^3 & )[x] == 0, x])]] > Chop[N[%]] > > -((2*(12 + 2^(1/3)*(27 - 3*Sqrt[69])^(2/3) + 6*3^(1/3)*(2/(9 - > Sqrt[69]))^(2/3))* > Log[1 + (2/(27 - 3*Sqrt[69]))^(1/3) + ((1/2)*(9 - > Sqrt[69]))^(1/3)/3^(2/3)])/ > (3*(6 + 2^(1/3)*(27 - 3*Sqrt[69])^(2/3) + 6*3^(1/3)*(2/(9 - > Sqrt[69]))^(2/3)))) - > (2*(8*I*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*2^(1/3)*3^(1/6)*(- > I + Sqrt[3]) + > (9 - Sqrt[69])^(1/3)*(-9 - 3*I*Sqrt[3] + I*Sqrt[23] + > Sqrt[69])))* > Log[1 + (-1 + I*Sqrt[3])/(2^(2/3)*(27 - 3*Sqrt[69])^(1/3)) - ((1 + > I*Sqrt[3])*((1/2)*(9 - Sqrt[69]))^(1/3))/(2*3^(2/3))])/ > (3*(4*I*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*2^(1/3)*3^(1/6)*(- > I + Sqrt[3]) + > (9 - Sqrt[69])^(1/3)*(-9 - 3*I*Sqrt[3] + I*Sqrt[23] + > Sqrt[69])))) - > (2*(8*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*I*2^(1/3)*3^(1/6)*(I > + Sqrt[3]) + > (9 - Sqrt[69])^(1/3)*(-9*I - 3*Sqrt[3] + Sqrt[23] + > I*Sqrt[69])))* > Log[1 + (-1 - I*Sqrt[3])/(2^(2/3)*(27 - 3*Sqrt[69])^(1/3)) + (I*(I > + Sqrt[3])*((1/2)*(9 - Sqrt[69]))^(1/3))/(2*3^(2/3))])/ > (3*(4*(9 - Sqrt[69])^(2/3) + > 2^(1/3)*3^(1/6)*(2*I*2^(1/3)*3^(1/6)*(I + Sqrt[3]) + > (9 - Sqrt[69])^(1/3)*(-9*I - 3*Sqrt[3] + Sqrt[23] + > I*Sqrt[69])))) > 0.3712169752602468 > > Am I missing something simpler than this setting? > Is it a built -in function that RootSum and this function acts like > the pair {RootReduce,ToRadicals} below? > > Reduce[1 + 2*x + 3*x^2 + x^3 == 0, Reals] > ToRadicals[%] > RootReduce[%] > > x == Root[1 + 2*#1 + 3*#1^2 + #1^3 & , 1] > x == -1 - (2/(3*(9 - Sqrt[69])))^(1/3) - ((1/2)*(9 - Sqrt[69]))^(1/3)/ > 3^(2/3) > x == Root[1 + 2*#1 + 3*#1^2 + #1^3 & , 1] > > > Thanks! > > Dimitris