Re: (Not trivial) Definite Integration of a rational function
- To: mathgroup at smc.vnet.net
- Subject: [mg74841] Re: (Not trivial) Definite Integration of a rational function
- From: "dimitris" <dimmechan at yahoo.com>
- Date: Sat, 7 Apr 2007 04:07:51 -0400 (EDT)
- References: <euslob$9nf$1@smc.vnet.net><ev502k$90n$1@smc.vnet.net>
> I think that when Mathematica integration algorithm deals with an
> improper integral
> does not use the Newton-Leibniz formula.
I think the following examples contribute to my statement...
Consider first the following integral
f = HoldForm[Integrate[1/(x^4 + 3*x^2 + 1)^8, {x, 0, Infinity}]]
Then
(*symbolic result*)
f // ReleaseHold // Timing
{48.843*Second, (21377637*Pi)/(160000000*Sqrt[5])}
(*numerical checking*)
{N[%[[2]]], ReleaseHold[f /. Integrate -> NIntegrate]}
{0.1877174241405516, 0.18771742413996437}
On the other hand
Integrate[1/(x^4 + 3*x^2 + 1)^8, x]
Timing[FullSimplify[Limit[%, x -> Infinity] - Limit[%, x -> 0]]]
(1/11200000000)*((1/(1 + 3*x^2 + x^4)^7)*(10*x*(921137791 +
10549280917*x^2 + 66382215753*x^4 + 260483882779*x^6 +
681502486770*x^8 + 1225622720238*x^10 + 1535886043461*x^12 +
1347185144039*x^14 + 826568154762*x^16 + 351855888230*x^18 +
101658554721*x^20 + 19007721747*x^22 + 2074781583*x^24 +
100424709*x^26)) +
63*Sqrt[2/(3 - Sqrt[5])]*(7970215 + 1530957*Sqrt[5])*ArcTan[Sqrt[2/
(3 - Sqrt[5])]*x] -
63*Sqrt[2/(3 + Sqrt[5])]*(-7970215 + 1530957*Sqrt[5])*ArcTan[Sqrt[2/
(3 + Sqrt[5])]*x])
{2.593999999999994*Second, (21377637*Pi)/(160000000*Sqrt[5])}
Plot[%%,{x,0,10}]
(*plot to be displayed*)
Of course someone could say that a considerable timing is spent in
order of check for singularities/convergence
and the stuff. But I don't think so
f /. Integrate[x___] :> Integrate[x, GenerateConditions -> False] //
ReleaseHold // Timing
{47.625*Second, (21377637*Pi)/(160000000*Sqrt[5])}
As another example consider
g = HoldForm[Integrate[Exp[-x] Log[x]^6, {x, 0, Infinity}]]
Then
Timing[ReleaseHold[g]]
{N[%[[2]]], ReleaseHold[g /. Integrate -> NIntegrate]}
{4.797000000000001*Second, EulerGamma^6 + (5*EulerGamma^4*Pi^2)/2 +
(9*EulerGamma^2*Pi^4)/4 + (61*Pi^6)/168 + 40*EulerGamma^3*Zeta[3] +
40*Zeta[3]^2 + 4*EulerGamma*(5*Pi^2*Zeta[3] + 36*Zeta[5])}
{715.0673625273188, 715.0673624721188}
On the other hand
Integrate[Exp[-x]*Log[x]^6, x]
Plot[%, {x, 0, 20}]
Timing[Limit[%%, x -> Infinity] - Limit[%%, x -> 0]]
x*(720*HypergeometricPFQ[{1, 1, 1, 1, 1, 1, 1}, {2, 2, 2, 2, 2, 2, 2},
-x] -
720*HypergeometricPFQ[{1, 1, 1, 1, 1, 1}, {2, 2, 2, 2, 2, 2}, -
x]*Log[x] +
360*HypergeometricPFQ[{1, 1, 1, 1, 1}, {2, 2, 2, 2, 2}, -
x]*Log[x]^2 - 120*HypergeometricPFQ[{1, 1, 1, 1}, {2, 2, 2, 2}, -x]*
Log[x]^3 + 30*HypergeometricPFQ[{1, 1, 1}, {2, 2, 2}, -x]*Log[x]^4 +
((-1 + E^x)*Log[x]^6)/(E^x*x) -
(6*Log[x]^5*(EulerGamma + Gamma[0, x] + Log[x]))/x)
(*plot to be displayed*)
{114.109*Second, Limit[x*(720*HypergeometricPFQ[{1, 1, 1, 1, 1, 1, 1},
{2, 2, 2, 2, 2, 2, 2}, -x] -
720*HypergeometricPFQ[{1, 1, 1, 1, 1, 1}, {2, 2, 2, 2, 2, 2}, -
x]*Log[x] +
360*HypergeometricPFQ[{1, 1, 1, 1, 1}, {2, 2, 2, 2, 2}, -
x]*Log[x]^2 -
120*HypergeometricPFQ[{1, 1, 1, 1}, {2, 2, 2, 2}, -x]*Log[x]^3 +
30*HypergeometricPFQ[{1, 1, 1}, {2, 2, 2}, -x]*Log[x]^4 + ((-1 +
E^x)*Log[x]^6)/(E^x*x) - (6*Log[x]^5*(EulerGamma + Gamma[0, x] +
Log[x]))/x), x -> Infinity]}
I=2Ee. the Newton-Leibniz formula failed to be applied.
Considering these examples (and others of course) and in view of the
relevant material appeared in the links I mentioned I hope I am right
and I don't fool myself!
Dimitris
=CF/=C7 dimitris =DD=E3=F1=E1=F8=E5:
> 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