Re: Discrepancy between Integrate and NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg83570] Re: Discrepancy between Integrate and NIntegrate
- From: "David W.Cantrell" <DWCantrell at sigmaxi.net>
- Date: Fri, 23 Nov 2007 05:31:49 -0500 (EST)
- References: <fhrr5o$552$1@smc.vnet.net>
chuck009 <dmilioto at comcast.com> wrote:
> I first use Integrate and then evaluate the results with n=3 and
> the limits t=pi and t=0. I then numerically evaluate the same
> integrand and the results are very different. Can anyone help
> me with this problem?
Sure. (BTW, I'm surprised that nobody has responded to your request yet.)
Your results don't _quite_ prove a bug. But I will show a bug (present in
version 5.2, and probably also in version 6) which is very closely related
to your discrepancy. I'll try to explain everything clearly.
> Clear[n]
> val1 = n*Pi*I*Integrate[(Cos[n*Pi*Exp[I*t]]/
> Cosh[n*Pi*Exp[I*t]])*Exp[I*t], t]
> N[(val1 /. {n -> 3, t -> Pi}) - (val1 /. {n -> 3, t -> 0})]
>
> n = 3;
> i1 = NIntegrate[n*Pi*I*(Cos[n*Pi*Exp[I*t]]/Cosh[n*Pi*Exp[I*t]])*
> Exp[I*t], {t, 0, Pi}]
>
> -1.2522017302861566 + 0.*I
> 7757.4213022137155 + 0.*I
As you might have supposed, the result of NIntegrate is correct. But let's
get in mind what it represents: the integral of Cos[z]/Cosh[z]
counterclockwise along the portion of the circle, centered at the origin
and having radius 3 Pi, in the upper half of the complex plane.
If we were to consider a closed curve C formed by the semicircle mentioned
above and its diameter along the real axis, then we can get the value of
the integral around C using the residue theorem. Within that contour, we
have only simple poles at I Pi/2, 3 I Pi/2, and 5 I Pi/2. Mathematica can
easily evaluate the residues. For example,
In[4]:= Residue[Cos[z]/Cosh[z], {z, 5 I Pi/2}]
Out[4]= (-I)*Cosh[(5*Pi)/2]
So then the value of the integral around C is
2 Pi I (-I Cosh[Pi/2] - (-I) Cosh[3*Pi/2] + (-I) Cosh[5*Pi/2])
which simplifies
In[7]:=
Simplify[2 Pi I (-I Cosh[Pi/2] - (-I) Cosh[3*Pi/2] + (-I)Cosh[5*Pi/2])]
Out[7]= 2*Pi*Cosh[Pi/2]*(1 - 2*Cosh[Pi])^2
In[8]:= N[%]
Out[8]= 7758.67
An important clue to your discrepancy is that the numerical value above is
the sum of the absolute values of your two numerical results, -1.25 and
7757.42 .
Now consider something which is closely related to your val1:
In[9]:=
N[Integrate[(Cos[3*Pi*E^(I*t)]/Cosh[3*Pi*E^(I*t)])*3*Pi*E^(I*t)*I,
{t, 0, Pi}]]
Out[9]= -1.2522 + 0.*I
This shows a bug, of a sort which is not uncommon in Mathematica. Now
you're probably wondering why I call that a bug, but had said earlier, when
you got that same numerical answer, that you hadn't quite proven a bug!
Well, it's because, although your val1 had discontinuites, you nonetheless
proceeded to evaluate the definite integral by using just F(b) - F(a), and
so it could be argued that you made the mistake, rather than Mathematica,
by not taking the discontinuites into account. But my Out[9] indicates that
Mathematica _itself_ makes the very same mistake, and so that is indeed a
bug.
We can, if careful, get a correct symbolic answer for the integral along
just your semicircle:
In[10]:=
2 Re[Integrate[(Cos[3*Pi*E^(I*t)]/Cosh[3*Pi*E^(I*t)])*3*Pi*E^(I*t)*I,
{t, Pi/2, Pi}]]
Out[10]=
2*Re[((-(1/4))*((2 + 2*I)*Hypergeometric2F1[1/2 - I/2, 1, 3/2 - I/2,
-E^(-6*Pi)] + (2 - 2*I)*Hypergeometric2F1[1/2 + I/2, 1, 3/2 + I/2,
-E^(-6*Pi)] + PolyGamma[0, 1/4 + I/4] + E^(6*Pi)*(PolyGamma[0, 1/4 - I/4] -
PolyGamma[0, 3/4 - I/4]) - PolyGamma[0, 3/4 + I/4]))/E^(3*Pi)]
In[11]:= N[%]
Out[11]= 7757.42
But I should explain what I did: I looked at graphs of real and imaginary
parts of Mathematica's indefinite integral. Although the real part had
discontinuities, it didn't have any between t = Pi/2 and Pi. Thus, by
symmetry, I knew that doubling the real part of the definite integral from
Pi/2 to Pi should correctly give the real part of the desired result.
Furthermore, by symmetry, the imaginary part of the integral from 0 to Pi/2
cancelled the part from Pi/2 to Pi. Thus, the real part, shown in Out[10],
is our complete symbolic answer.
Finally, note that Mathematica correctly calculates the integral along the
real axis from -3 Pi to 3 Pi:
In[15]:= Integrate[Cos[x]/Cosh[x], {x, -3 Pi, 3 Pi}]
Out[15]=
((1/2)*I*(Beta[-E^(-6*Pi), 1/2 - I/2, 0] - Beta[-E^(6*Pi), 1/2 - I/2, 0] +
E^Pi*(Beta[-E^(-6*Pi), 1/2 + I/2, 0] - Beta[-E^(6*Pi), 1/2 + I/2,
0])))/E^(Pi/2)
In[16]:= N[%]
Out[16]= 1.2522 + 0.*I
Then adding Out[10] and Out[15], we get a result in numerical agreement
with Out[7], our integral around C which we had calculated using the
residue theorem.
David W. Cantrell