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