[Date Index]
[Thread Index]
[Author Index]
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
Prev by Date:
**Re: How to make results from Integrate and NIntegrate agree**
Next by Date:
**Differential equations**
Previous by thread:
**How to make results from Integrate and NIntegrate agree**
Next by thread:
**Re: How to make results from Integrate and NIntegrate agree**
| |