Re: How to make results from Integrate and NIntegrate agree

*To*: mathgroup at smc.vnet.net*Subject*: [mg66970] Re: How to make results from Integrate and NIntegrate agree*From*: ab_def at prontomail.com*Date*: Mon, 5 Jun 2006 03:48:46 -0400 (EDT)*References*: <200605240701.DAA08462@smc.vnet.net><e53leb$3t1$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Andrzej Kozlowski wrote: > 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 In fact, Integrate can handle integrands that can be expanded into a finite number of cases (not greater than $MaxPiecewiseCases) as well as some simple infinite-piecewise functions (e.g., some integrals involving UnitStep[Sin[t]]). For example: In[1]:= PiecewiseExpand[FractionalPart[1/t], t > 1/2] Out[1]= Piecewise[{{1/t, t > 1}}, (1 - t)/t] In[2]:= Assuming[x > 1/2, Integrate[FractionalPart[1/t], {t, 1/2, x}]] Out[2]= Piecewise[{{(1/2)*(-1 + 2*Log[2] + 2*Log[x]), x > 1}}, (1/2)*(1 - 2*x + 2*Log[2*x])] PiecewiseIntegrate doesn't know about the singularities of ArcTan, it will simply pass it on to Integrate. To get PiecewiseIntegrate to find the antiderivative of FractionalPart[1/x] on the whole real line, we need to deal with more problems in sums of logarithms: In[3]:= Sum[Log[(k + 2)/(k + 1)] - 1/(k + 2), {k, 0, Infinity}] Out[3]= Indeterminate In[4]:= Sum[Log[(k + 2)/(k + 1)] - 1/(k + 2), {k, a, b}] Out[4]= Indeterminate So first we add some rules for such sums: Unprotect[Sum]; logF = True; Sum[r_.*Log[e_], {k_, a_, b_}] := Module[{p, q}, {p, q} = Coefficient[e, k, {1, 0}]; r*Log[p^(b + 1 - a)*Gamma[b + 1 + q/p]/Gamma[a + q/p]] ] /; (!NumericQ[a] || !NumericQ[b]) && PolynomialQ[e, k] && Exponent[e, k] == 1 && FreeQ[r, k] Sum[e_, {k_, a_, Infinity | DirectedInfinity[1]}] := Module[{n}, Limit[Sum[e, {k, a, n}], n -> Infinity]] /; !FreeQ[e, Log] Sum[e_, {k_, a_, b_}] := Block[{logF = False}, If[Head@ # === Plus, Sum[#, {k, a, b}]& /@ #, Sum[#, {k, a, b}] ]&@ Apart[e, k] ] /; (!NumericQ[a] || !NumericQ[b]) && logF && !FreeQ[e, Log] && !MemberQ[{a, b}, _DirectedInfinity] And finally In[10]:= PiecewiseIntegrate[FractionalPart[1/t], {t, 0, x}] // PiecewiseExpand // Simplify Out[10]= Piecewise[{{1 - EulerGamma, x == -1 || x == 1}, {2 - EulerGamma + x + Log[-x], -1 < x < -(1/2)}, {1 - EulerGamma + Log[-x], x < -1}, {1 - EulerGamma + Log[x], x > 1}, {2 - EulerGamma - x + Log[x], 1/2 < x < 1}, {(-(1/2))*(-2 + x)*Ceiling[1/x] - (1/2)*x*Ceiling[1/x]^2 - 1/Floor[1/x] + (1/2)*(-2 + x)*Floor[1/x] + (1/2)*x*Floor[1/x]^2 - Log[Gamma[1 - Floor[1/x]]] + Log[((-x)^(Ceiling[1/x] - Floor[1/x])*Gamma[1 - Floor[1/x]])/Gamma[1 - Ceiling[1/x]]] + Log[Gamma[-Floor[1/x]]] + PolyGamma[0, 1 - Ceiling[1/x]] - PolyGamma[0, 1 - Floor[1/x]] + PolyGamma[0, -Floor[1/x]], Inequality[-(1/2), LessEqual, x, Less, 0]}}, 1 + 1/Floor[1/x] - x*Floor[1/x] + Log[Gamma[1 + Floor[1/x]]] - Log[Gamma[2 + Floor[1/x]]] + Log[(x*Gamma[2 + Floor[1/x]])/Gamma[1 + Floor[1/x]]] + PolyGamma[0, Floor[1/x]]] The above rule for the sum of Log[a*k + b] is sufficient here but it's not valid generically. It's often the case that closed forms for sums are not valid for all values of the parameters: In[11]:= Sum[(-k)^(1/2), {k, a, b}] Out[11]= I*(Zeta[-(1/2), a] - Zeta[-(1/2), 1 + b]) (Not correct for a < b < 0). Also there are some unrelated problems with piecewise functions, sometimes rather strange ones: In[12]:= PiecewiseExpand[Floor[1/Sqrt[x]], 0 < x < 1/Sqrt[2]] Out[12]= 0 In[13]:= Limit[Floor[1/x], x -> 0] Out[13]= 100000000 So it's always a good idea to double-check the piecewise results. Maxim Rytin m.r at inbox.ru