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