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