Re: How to make results from Integrate and NIntegrate agree
- To: mathgroup at smc.vnet.net
- Subject: [mg66650] Re: How to make results from Integrate and NIntegrate agree
- From: "David W.Cantrell" <DWCantrell at sigmaxi.org>
- Date: Thu, 25 May 2006 02:58:56 -0400 (EDT)
- References: <e510la$8di$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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. I'm sorry that I can't answer your questions. Nonetheless, I hope you find my comments helpful. I seems to me that what you need is a _continuous_ antiderivative. Your question here inspired me to write "A continuous antiderivative for the floor function", which I posted to sci.math a little while ago. Please see <http://groups.google.com/group/sci.math/msg/bfd05d9f92fff65d>. I think you'll find it to be pertinent. Of course, your specific function FrPart[1/x] is messier than Floor[x]. Nonetheless, it happens that, for FrPart[1/x], an antiderivative which is continuous for x > 0 can be expressed in closed form in Mathematica: Mod[1, x] + PolyGamma[Floor[1/x] + 1] + Log[x] Unfortunately, I don't know how to force Mathematica to give such an antiderivative. David W. Cantrell