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: [mg66694] Re: How to make results from Integrate and NIntegrate agree
  • From: Maxim <m.r at inbox.ru>
  • Date: Sat, 27 May 2006 03:52:08 -0400 (EDT)
  • References: <e510la$8di$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On Wed, 24 May 2006 07:06:18 +0000 (UTC), <riazzi at hotmail.com> wrote:

> Hello,
>
> I need to evaluate the integral of a discontinuous function for which
> the results from Integrate and NIntegrate do not agree.   A more simple
> function which demonstrates the problem I face is the following, which
> is a simplified form (for nonnegative, real arguments) of the
> FractionalPart function from functions.wolfram.com.  I'm using this
> form because Integrate returns my integral unevaluated when I use the
> built in FractionalPart[].
>
> In[1]=
> FrPart[x_]:=(1/2)-ArcTan[Cot[Pi*x]]/Pi
>
>
> So here's the example function to be integrated, showing it's
> discontinuities at {1/n : n a Natural number}:
>
> In[2]=
> Plot[FrPart[1/x], {x,1/6,1},PlotPoints->100];
>
>
> If I integrate this numerically, I get what appears to be the right
> answer, and the plot looks like what I would expect to see:
>
> In[3]=
> nintgrl[x_] := If[x = (1/6), 0, NIntegrate[FrPart[1/u], {u, (1/6),
> x},
>     MaxRecursion -> 30, WorkingPrecision -> 30, PrecisionGoal -> 25]]
> Plot[nintgrl[x], {x, 1/6, 1}, PlotStyle -> {Red}];
> nintgrl[1]
>
> Out[3]= <<Plot omitted>>
> 0.3417594692
>
> If I evaluate the integral piecewise, I obtain an answer that agrees
> with NIntegrate's:
>
> In[4]=
> area = Sum[Integrate[1/u -n, {u,1/(n+1),1/n}],{n,1,5}]
> N[area]
>
> Out[4]= -(29/20) + Log[6/5] + Log[5/4] + Log[4/3] + Log[3/2] + Log[2]
> 0.341759
>
>
> If I evaluate the integral with Integrate, here is what I get:
>
> In[5]=
> FullSimplify[Integrate[FrPart[1/x], x]]
>
> Out[5]= (x/2) - (x * ArcTan[Cot[Pi/x]])/Pi + Log[x] - 1
>
> (NOTE: I am aware that having Integrate evaluate the definite integral
> Integrate[FrPart[1/u], {u,1/6,1}] returns the correct answer in this
> case.  For my more complicated function over a different range it does
> not.)
>
>
> Here's how the Integrate result looks, with a constant of integration
> of 5/6 + Log[6] added to make the plot start at zero like the
> NIntegrate result:
>
> In[6]=
> Plot[(x/2) - (x * ArcTan[Cot[Pi/x]])/Pi + Log[x] - 1 + 5/6 + Log[6],
> {x, 1/6, 1}, PlotStyle -> Blue];
>
>
> Looking at limits as x approaches (1/n) from either direction, I see
> that the magnitude of the jump at (1/n) equals (1/n).  If I remove the
> jumps (by subtracting off their contribution via the function adj[x]),
> the NIntegrate and Integrate plots overlap perfectly:
>
> In[7]=
> adj[x_] := Sum[(1/n)*UnitStep[x - 1/n], {n, 2,  5}];
> adjPlot = Plot[(x/2) - (x * ArcTan[Cot[Pi/x]])/Pi + Log[x] - 1 + 5/6 +
> Log[6] - adj[x], {x, 1/6, 1}, PlotStyle -> Blue];
>
>
> But I don't want to subtract off the contribution of the jumps, nor do
> I want to sum up a (for me infinite) collection of piecewise terms.  So
> what I'm looking for is a way to make Integrate produce a closed-form
> expression for my integrand that does not produce these
> discontinuities, plus some insight as to why they appear.  I'd guess
> the latter has to do with branch cuts in the complex plane for the
> functions involved, and where the solution might have to do with
> picking a different contour or somesuch, but that's getting into areas
> I've long forgotten or never learned.
>
> Thanks in advance,
>     Will R.
>

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).

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

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

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: parenthesis in reactions...
  • Next by Date: Re: Newbie: Rotataion2D use problem
  • Previous by thread: Re: How to make results from Integrate and NIntegrate agree
  • Next by thread: Re: How to make results from Integrate and NIntegrate agree