Services & ResourcesWolfram ForumsMathGroup Archive

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

How to make results from Integrate and NIntegrate agree

• To: mathgroup at smc.vnet.net
• Subject: [mg66610] How to make results from Integrate and NIntegrate agree
• From: riazzi at hotmail.com
• Date: Wed, 24 May 2006 03:01:44 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

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.

• Follow-Ups:
  • Re: How to make results from Integrate and NIntegrate agree
    • From: Andrzej Kozlowski <akoz@mimuw.edu.pl>