Re: (Not trivial) Definite Integration of a rational function
- To: mathgroup at smc.vnet.net
- Subject: [mg74755] Re: (Not trivial) Definite Integration of a rational function
- From: "Michael Weyrauch" <michael.weyrauch at gmx.de>
- Date: Wed, 4 Apr 2007 04:00:35 -0400 (EDT)
- References: <euslob$9nf$1@smc.vnet.net>
Hello, well, another one of your mystery stories involving Integrate. It is obvious that the only "problem" to evaluate this integral is the calculation of the limit of the antiderivative for x going to Infinity. After some serious thinking (more than a minute on my PC) Mathematica 5.2 returns 0 for this limit, which is correct. (Mathematica 4.2 returns Infinity which is wrong.) The limit is calculated from the sum of a number of logarithms which separately all diverge for x-> Infinity, of course. Explicitly the antiderivative reads 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) & ]), and if I was to calculate the limit for x-> Infinity I would just neglect the various constants under the Logs and find 2*(-Log[x] + RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[x] + 2*Log[x]*#1 + Log[x]*#1^2)/ (2 + 6*#1 + 3*#1^2) & ]) which Mathematica simplifies to 0 for all x immediately. So, the calculation of the limit is actually a matter of milliseconds (for a CAS), and why Mathematica does not really manage to do this in practise without the explicit "help" remains a mystery. It more and more appears to me that one of the issues with Integrate is the calculation of limits. As we have already seen in other examples, it is this step which for obscure reasons (to me) seems to be particularly difficult. The questions to be answered are: What are the algorithms and methods Mathematica uses for Limit[], and why is it unable to use simple steps (tricks) like the one used above. A real mystery, of course, is the hangup Mathematica runs into, if we ask to calculate the definite integral directly. I guess, this again can only be answered by the developers... Michael Weyrauch "dimitris" <dimmechan at yahoo.com> schrieb im Newsbeitrag news:euslob$9nf$1 at smc.vnet.net... >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 > >