Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: integrate

  • To: mathgroup at smc.vnet.net
  • Subject: [mg72453] Re: integrate
  • From: "David W. Cantrell" <DWCantrell at sigmaxi.net>
  • Date: Thu, 28 Dec 2006 05:36:00 -0500 (EST)
  • References: <emgg3i$2ca$1@smc.vnet.net>

"dimitris" <dimmechan at yahoo.com> wrote:
For your third question:
> Integrate[Sqrt[x^2 + 1]/Sqrt[1 + x^6], {x, 0, 10}]
> (-(-1)^(1/6))*EllipticF[I*ArcSinh[10*(-1)^(1/3)], (-1)^(2/3)]
>
> N[%]
> -0.10016641038463325 - 1.6857503548125956*I
>
> NIntegrate[Sqrt[x^2 + 1]/Sqrt[1 + x^6], {x, 0, 10}]
> 2.056349237110889
>
> --------> Any ideas to "help" Mathematica to give a correct answer?

I have an answer, but getting there was not pleasant. I hope others
will have better suggestions.

Part 1:
I wanted to know why the symbolic answer was wrong, so I looked at.
In[7]:= Integrate[Sqrt[(1 + x^2)/(1 + x^6)], x]
Out[7]=
(-(-1)^(1/6))*Sqrt[1 - (-1)^(1/3)*x^2]*Sqrt[1 + (-1)^(2/3)*x^2]*
Sqrt[1/(1 - x^2 + x^4)]*EllipticF[I*ArcSinh[(-1)^(1/3)*x], (-1)^(2/3)]

The product of the algebraic factors involving x should be just 1, but
I was not able to get Mathematica to give that:

In[8]:=
FullSimplify[Sqrt[1 - (-1)^(1/3)*x^2]*Sqrt[1 + (-1)^(2/3)*x^2]*
Sqrt[1/(1 - x^2 + x^4)], x > 0]
Out[8]=
Sqrt[1 - (-1)^(1/3)*x^2]*Sqrt[(1 + (-1)^(2/3)*x^2)/(1 - x^2 + x^4)]

Maybe I didn't try hard enough. Can someone suggest how I should get
the above to simplify to 1, perhaps assuming x > 0? Anyway, I then
simplified Out[7] myself to

In[9]:= -(-1)^(1/6)*EllipticF[I*ArcSinh[(-1)^(1/3)*x], (-1)^(2/3)]

and then plotted the real part of that from x = 0 to 10. A jump
discontinuity is apparent, at roughly x = 2.45. If we ignore that
discontinuity and try to evaluate the integral from x = 0 to 10 by the
Fundamental Theorem, we get the same wrong answer that Mathematica
got. So now I think I know what Mathematica did wrong. But I'm a bit
surprised. Even when Mathematica returns antiderivatives having such
discontinuities, it normally takes such discontinuities correctly into
account when doing _definite_ integrals. But it didn't do so here, so
there is indeed a bug to be exterminated.

Part 2:
I thought I would help Mathematica by taking the discontinuity, caused
by a branch cut, into account by hand. But unfortunately it wasn't
obvious exactly where the cut was. So I went to the Wolfram Functions
site and looked under EllipticF. It was not helpful. The only relevant
statement was "Branch cut locations: complicated." I then gave up on
determining the exact location. But if I could have done that, it
would have been nicer than what follows. (Anyone know exactly where
the discontinuity is?)

Part 3:
I sought to transform the original integral into another one which
Mathematica might handle better. My first thought was to use the
substitution u = x^2, which led to

In[16]:= 1/2*Integrate[Sqrt[(1 + u)/(u*(1 + u^3))], u]
Out[16]=
(-1)^(1/6)*Sqrt[1 - (-1)^(1/3)/u]*Sqrt[1 + (-1)^(2/3)/u]*u^(3/2)*
Sqrt[1/(u - u^2 + u^3)]*EllipticF[I*ArcSinh[(-1)^(1/3)/Sqrt[u]],(-1)^(2/3)]

Simplifying by hand and putting the result back in terms of x gave

In[17]:= (-1)^(1/6)*EllipticF[I*ArcSinh[(-1)^(1/3)/x], (-1)^(2/3)].

Although different from In[9], this result still has a discontinuity.
But thankfully, it's in a different place, roughly x = 0.4. That
allows us to combine In[9] and In[17], using the former to integrate
from x0 to 1 and the latter to integrate from 1 to x1. Thus, assuming
x0 < 1 < x1, we can state that

Integrate[Sqrt[(1 + x^2)/(1 + x^6)], {x, x0, x1}]

is

In[21]:=
int = (-1)^(1/6)*(-2*EllipticF[I*ArcSinh[(-1)^(1/3)], (-1)^(2/3)] +
EllipticF[I*ArcSinh[(-1)^(1/3)*x0], (-1)^(2/3)] +
EllipticF[I*ArcSinh[(-1)^(1/3)/x1], (-1)^(2/3)]);

In[22]:= int /. {x0 -> 0, x1 -> 10}
Out[22]=
(-1)^(1/6)*(EllipticF[I*ArcSinh[(1/10)*(-1)^(1/3)], (-1)^(2/3)] -
2*EllipticF[I*ArcSinh[(-1)^(1/3)], (-1)^(2/3)])

which is the precise symbolic result desired.

In[23]:= N[%]
Out[23]= 2.0563492371150103 + 3.3306690738754696*^-16*I

David W. Cantrell


  • Prev by Date: Re: Please help overcome problems with integration
  • Next by Date: Re: Any simple way to flatten all but the bottom level?
  • Previous by thread: integrate
  • Next by thread: Re: Re: integrate