Re: Integrate vs NIntegrate

*To*: mathgroup at smc.vnet.net*Subject*: [mg89145] Re: Integrate vs NIntegrate*From*: "David W.Cantrell" <DWCantrell at sigmaxi.net>*Date*: Tue, 27 May 2008 07:13:42 -0400 (EDT)*References*: <g1e35c$er3$1@smc.vnet.net>

"Armen Kocharyan" <armen.kocharyan at gmail.com> wrote: > Dear group, > > Why do I get different results (Res1, Res2) in Mathematica 6.0.1.0? > > *q=-1/2; a=0; b=3;* > > *h[x_]=(1+x^3)^q;* > > *f[x_]=Integrate[h[x],x];* > > *Res1=N[f[b]-f[a]]* Res1 will not be the desired value. Reason: f is not an antiderivative of h on the interval [0, 3] due to the fact that f has a discontinuity at x = 2. But, knowing that f has a discontinuity there, we should be able to use f to calculate the desired value. Namely, the desired value should be given by f[b] - Limit[f[x], x -> 2, Direction -> -1] + Limit[f[x], x -> 2, Direction -> 1] - f[a] Strangely, that doesn't work either! The reason that the above does not give the desired value is that Mathematica has a _bug_ which make it think that the first limit is the same as the second limit. We can, nonetheless, use f to get a reasonably accurate value for the integral: In[20]:= N[f[b] - f[2 + 10^-7] + f[2 - 10^-7] - f[a]] Out[20]= 1.65267 + 0. I BTW, although Mathematica's f, which is in terms of EllipticF, has a discontinuity at x = 2, one can work out an antiderivative in terms of EllipticF which is valid on [0, 3]: In[21]:= fcont[x_] = 1/3^(1/4)* EllipticF[ArcCos[(Sqrt[3] - 1 - x)/(Sqrt[3] + 1 + x)], (2 + Sqrt[3])/4]; Of course, that allows us to calculate the definite integral more simply, in the way you had intended: In[22]:= N[fcont[b] - fcont[a]] Out[22]= 1.65267 > *Res2=NIntegrate[h[x],{x,a,b}]* > > *Res1 = -2.55387-2.42865 i* > > *Res2 = 1.65267* > > I'm assuming that Res2 is the correct answer. Right. Another method: If you prefer to have things in terms of Hypergeometric2F1, instead of EllipticF, then In[23]:= Assuming[xf > 0, Integrate[(1 + x^3)^(-1/2), {x, 0, xf}]] Out[23]= xf*Hypergeometric2F1[1/3, 1/2, 4/3, -xf^3] In[24]:= N[% /. xf -> 3] Out[24]= 1.65267 David W. Cantrell