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