MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: How to make results from Integrate and NIntegrate agree

  • To: mathgroup at smc.vnet.net
  • Subject: [mg66728] Re: How to make results from Integrate and NIntegrate agree
  • From: "David W.Cantrell" <DWCantrell at sigmaxi.org>
  • Date: Sat, 27 May 2006 21:04:02 -0400 (EDT)
  • References: <e510la$8di$1@smc.vnet.net> <e592lv$3c0$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Maxim <m.r at inbox.ru> wrote:
> PiecewiseIntegrate
> ( http://library.wolfram.com/infocenter/MathSource/5117/ ) can find a
> continuous antiderivative of FractionalPart[x] on the real line:
>
> In[2]:= PiecewiseIntegrate[FractionalPart[t], {t, 0, x}]
>
> Out[2]= If[x == -2, 1/2, 0] + If[x == 2, 1/2, 0] + If[-2 < x < -1,
> (1/2)*(1 + 2*x + x^2), 0] + If[-1 < x < 0, x^2/2, 0] + If[0 < x < 1,
> x^2/2, 0] + If[1 < x < 2, (1/2)*(1 - 2*x + x^2), 0] + If[2 < x,
> (1/12)*(-6 + Ceiling[x] + 6*x*Ceiling[x] + 6*x^2*Ceiling[x] -
> 3*Ceiling[x]^2 - 6*x*Ceiling[x]^2 + 2*Ceiling[x]^3 + 5*Floor[x] -
> 6*x*Floor[x] - 6*x^2*Floor[x] + 3*Floor[x]^2 + 6*x*Floor[x]^2 -
> 2*Floor[x]^3), 0] + If[x < -2, (1/12)*(-6 + Ceiling[-x] - 6*x*Ceiling[-x]
> + 6*x^2*Ceiling[-x] - 3*Ceiling[-x]^2 + 6*x*Ceiling[-x]^2 +
> 2*Ceiling[-x]^3 + 5*Floor[-x] + 6*x*Floor[-x] - 6*x^2*Floor[-x] +
> 3*Floor[-x]^2 - 6*x*Floor[-x]^2 - 2*Floor[-x]^3), 0] + If[1 <= x, 1/2, 0]
> + If[x <= -1, 1/2, 0]
>
> http://functions.wolfram.com/04.05.21.0005.01 gives a simpler
> antiderivative which is continuous on (-1, Infinity).

Right, their evaluation of the definite integral is not valid for a < -1.
Their entry should be modified either to indicate that restriction or,
better, to give a result which is valid for all real a.

Earlier in this thread, I posted a link to a thread in which I gave a
continuous antiderivative for Floor[x] on R. Since FractionalPart[x] is
easily written in terms of Floor[x], we can easily get a neat result for
Integrate[FractionalPart[t], {t, 0, a}] which is valid for all real a:

(FractionalPart[a]^2 - Sign[a]*FractionalPart[a] + Abs[a])/2

> Integration of FractionalPart[1/x] runs into some bugs in sums of
> logarithms:

Earlier in this thread, I gave a continuous antiderivative for that, valid
for x > 0. Of course, that result can be modified trivially to get an
antiderivative which is continuous on R,
namely, Integrate[FractionalPart[1/t], {t, 0, x}] is

If[x==0, 0, Mod[1, Abs[x]] + PolyGamma[Floor[1/Abs[x]] + 1] + Log[Abs[x]]]

> In[3]:= Sum[Log[2*k] - Log[2], {k, a, b}]
>
> Out[3]= (-(1 - a + b))*Log[2] - Log[Gamma[a]] + Log[Gamma[1 + b]]
>
> The problem arises when the summand includes terms of the form Log[p*k]
> (Mathematica loses the factor p), so we add a rule to expand the
> logarithms of products:
>
> In[4]:= Unprotect[Log];
>    Log[a_*x_] := Log[x] + Log[a] /; Refine[Positive[a]]
>
> In[6]:= int = PiecewiseIntegrate[FractionalPart[1/t], {t, eps, 1},
>        Assumptions -> 0 < eps < 1];
>
> Then the integral comes out right but a sum is left unevaluated; we can
> find the closed form if we perform the summation term by term:
>
> In[7]:= int = int /. s : HoldPattern[Sum[e_, {k_, a_ : 1, b_}]] :>
>      (If[Head@ # =!= Plus, s, Sum[#, {k, a, b}]& /@ #]&@ Apart[e, k]) //
>    PiecewiseExpand
>
> Out[7]= Piecewise[{{(-1 + 2*Log[2])/2, eps == 1/2}, {(-5 + 6*Log[3/2]
> + 6*Log[2])/6, eps == 1/3}, {(-3 + 4*eps - 2*Log[eps])/2, Inequality[1/3,
> Less, eps, Less, 1/2]}, {-1 + eps - Log[eps], eps > 1/2}}, (2 -
> 2*EulerGamma - 2*Ceiling[eps^(-1)] - eps*Ceiling[eps^(-1)]
> + eps*Ceiling[eps^(-1)]^2 + 2*Floor[eps^(-1)] + eps*Floor[eps^(-1)] -
> eps*Floor[eps^(-1)]^2 - 2*Ceiling[eps^(-1)]*Log[eps]
> + 2*Floor[eps^(-1)]*Log[eps] - 2*Log[Gamma[Ceiling[eps^(-1)]]]
> + 2*Log[Gamma[1 + Floor[eps^(-1)]]] - 2*PolyGamma[0, 1
> + Floor[eps^(-1)]])/2]
>
> In[8]:= int /. eps -> 1/6 // Simplify
>
> Out[8]= -(29/20) + Log[6]
>
> We can also find the limit at 0; when the limit exists it will coincide
> with the limit over the sequence eps = 1, 1/2, 1/3, ...:
>
> In[9]:= Limit[Refine[int, Element[1/eps, Integers]], eps -> 0]
>
> Out[9]= 1 - EulerGamma

Right. Or more simply, just substitute 1 for x in the definite integral I
gave above. That gives just PolyGamma[2], which is indeed 1 - EulerGamma.

Regards,
David W. Cantrell


  • Prev by Date: RE: Different text style for labels and ticks
  • Next by Date: Re: Interval[{a,b}]-Interval[{a,b}] = 0?
  • Previous by thread: Re: How to make results from Integrate and NIntegrate agree
  • Next by thread: Prefix function syntax f @ x