Re: How to make results from Integrate and NIntegrate agree

*To*: mathgroup at smc.vnet.net*Subject*: [mg66642] Re: [mg66610] How to make results from Integrate and NIntegrate agree*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Thu, 25 May 2006 02:58:21 -0400 (EDT)*References*: <200605240701.DAA08462@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

On 24 May 2006, at 16:01, 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. > I suggest downloading from MathSource Maxim Rytin's excellent package "Piecewise". http://library.wolfram.com/infocenter/MathSource/5117/ It can deal correctly with FractionalPart as well as with your function FrPart: << Piecewise` Simplify[PiecewiseIntegrate[FractionalPart[1/x], {x, 1/6, 1}]] -(29/20) + Log[6] N[%] 0.341759 FrPart[x_] := 1/2 - ArcTan[Cot[Pi*x]]/Pi PiecewiseIntegrate[FrPart[1/x], {x, 1/6, 1}] -(29/20) + Log[6] The package also contains a numerical integrator NPiecewiseIntegrate: NPiecewiseIntegrate[FractionalPart[1/x], {x, 1/6, 1}] 0.3417594692280566 I am not sure if PiecewiseIntegrate can manage a general formula for PiecewiseIntegrate[FractionalPart[1/x], {x, a, b}] It seems to me a lot to ask and as it appears to be taking a long time, which I could not to spare today I just aborted. Andrzej Kozlowski Tokyo, Japan

**References**:**How to make results from Integrate and NIntegrate agree***From:*riazzi@hotmail.com