       Re: [Mathematica 6] Integrate strange result

• To: mathgroup at smc.vnet.net
• Subject: [mg78376] Re: [mg78325] [Mathematica 6] Integrate strange result
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Fri, 29 Jun 2007 05:51:18 -0400 (EDT)
• References: <200706280829.EAA20660@smc.vnet.net> <BBF54CC9-1358-4764-BA31-BD8856F70C79@mimuw.edu.pl>

```I wrote the staff below very late at night, which not surprisingly,
proved a foolish thing to do. Had I gone to bed instead I would have
been better rested today and would not have posted a lot of nonsense
that I am having to correct now.

First of all the last bit - I simply wrote EllipticK instead of
EllipticE and "discovered" a bug in Mathematica. This is more or less
the same kind of discovery as used to be sometimes made by amateur
scientists who neglected to clean their microscope (or telescope).
Anyway:

N[2*(EllipticE[-k] + Sqrt[k + 1]*EllipticE[
k/(k + 1)]) /. k -> 5, 10]
11.32079298538350925`10.000000000000004

N[4*EllipticE[-k] /. k -> 5, 10]
11.32079298538350925`10.000000000000004

So everything is fine.

Let's return to the original problem and try to analyze it a bit:

EllipticE[m] is by definition the integral:

Integrate[Sqrt[1 - m*Sin[t]^2], {t, 0, Pi/2}]

and has a branch discontinuity in the complex plane running from 1
to Infinity. We have:

Integrate[Sqrt[1 + m*Sin[t]^2], {t, 0, 2Pi}] == Integrate[Sqrt[1 +
m*Sin[t]^2], {t, 0, Pi/2}]+Integrate[Sqrt[1 + m*Sin[t]^2], {t, Pi/2,
Pi}]+Integrate[Sqrt[1 + m*Sin[t]^2], {t, Pi, 3Pi/2}]+Integrate[Sqrt[1
+ m*Sin[t]^2], {t, 3Pi/2, 2Pi}]

It is immediately obvious (look at the graph of Sin[t]^2) that all
theses integrals are equal and all are equal to EllipticE[-m], at
least when m does not lie on the branch cut, in other words we do not
have Im[m]==0 && Re[-m] >=1, that is Im[m]!=0 || Re[m]>-1. So the answer

Integrate[f, {t, 0, 2*Pi}]

If[Re[k] > -1, 2*(EllipticE[-k] + Sqrt[k + 1]*EllipticE[k/(k + 1)]),
Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0, 2*Pi},
Assumptions -> Re[k] <= -1]]

is correct, but is only a "part of the truth". I think (unless I am
again confused) a better answer would be:

If[Re[k] > -1 || Im[k] !=0 ,4*EllipticE[-k],
Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0, 2*Pi},
Assumptions -> Re[k] <= -1]]

In fact 4*EllipticE[-k] == 2*(EllipticE[-k] + Sqrt[k + 1]*EllipticE
[k/(k + 1)]) at least outside the branch cut, but FullSimplify seems
unable to prove it. We conclude that the answer Mathematica returns to

Assuming[Element[k, Reals], Integrate[f, {t, 0, 2*Pi}]]

If[k >= -1, 4*EllipticE[-k], Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0,
2*Pi}, Assumptions -> k < -1]]

is perfectly correct. So the original question reduces to "what about
integrating along the branch cut of EllipticE", that is what about:

f = (1 + k*Sin[t]^2)^(1/2);
Assuming[Element[k, Reals] && k < -1,   Integrate[f, {t, 0,
2*Pi}]]

?

The answer 0 given by Mathematica 6, is clearly wrong under any
interpretation and must be considered a bug. But should Mathematica
give any answer in this situation? Well, perhaps, because  for a
particular case where m lies on the branch cut of EllipticE
Mathematica returns:

Integrate[(1 - 2*Sin[t]^2)^(1/2), {t, 0, 2*Pi}]

4*EllipticE

It seems to me that with using a definition of EllipticE along the
branch cut given by semi-continuity this answer is correct, and so is
4*EllipticE[-k] for the integral without any restrictions on k (I
perhaps that is the answer Integrate should return (it does so if you
Integrate with a general limit b and then use the replacement rule b -
> 2Pi).

Anyway, my conclusion is: we have found one bug is one bug in
Integrate (the 0 answer in the original post) and the possibility
that the answer to (Integrate[f, {t, 0, 2*Pi}]) could be more general
but other than that everything seems fine.

Andrzej Kozlowski

On 28 Jun 2007, at 23:14, Andrzej Kozlowski wrote:

> *This message was transferred with a trial version of CommuniGate
> (tm) Pro*
>
> On 28 Jun 2007, at 17:29, Nasser Abbasi wrote:
>
>> In:= \$Version
>> Out= 6.0 for Microsoft Windows (32-bit) (April 27, 2007)
>>
>> The problem is this: When I ask Mathematica to integrate something
>> involving some constant parameter, it gives an answer for some range
>> of the constant involved, say range A, and it says there is no answer
>> for the other range, say range B
>>
>> But next, when I specify, using assumptions, that the constant is now
>> in range B, then Mathematica does given an answer to the integral !
>>
>> Why does it in one case say there is no answer if the constant is in
>> range B, and in the second case it gives an answer when the constant
>> is in range B?
>>
>> case 1
>> ========
>>
>> In:= f = (1 + k*Sin[a]^2)^(1/2);
>> In:= Assuming[Element[k, Reals], Integrate[f, {a, 0, 2*Pi}]]
>>
>> Out= If[k >= -1, 4*EllipticE[-k],   Integrate[Sqrt[1 +
>> k*Sin[a]^2],    {a, 0, 2*Pi}, Assumptions ->     k < -1]]
>>
>> notice in the above, it says for k<-1 there is NO answer
>>
>> case 2
>> =======
>> In:= f = (1 + k*Sin[a]^2)^(1/2);
>>>> In:= Assuming[Element[k, Reals] && k < -1,   Integrate[f,
>>>> {a, 0,
>>>> 2*Pi}]]
>>>>
>>>> Out= 0
>>
>> Notice in the above, for k<-1  it gives an answer, which is zero.
>>
>> What is it I am missing here?
>>
>> thanks
>> Nasser
>>
>>
>
>
> The first thing you are missing here is that the second answer is
> wrong:
>
> Integrate[(1 - 2*Sin[a]^2)^(1/2), {a, 0, 2*Pi}]
>
> 4*EllipticE
>
> This is, as one would expect, a non-real number:
>
> N[%, 5]
> 2.3962804693361682454`5.000000000000002 +
>   2.3962804693361682454`5.000000000000002*I
>
> When Mathematica returns an answer in this conditional form it is
> not actually saying that when the condition is not satisfied "there
> is no answer". It is only saying that it could not establish an
> that is valid under the given conditions, but this does not mean
> that the answer (or some other answer) is not valid when this
> condition is not satisfied. All it means that Integrate has
> returned a set f  conditions under which it "believes" it has been
> able to establish that the answer it gave is valid.
>
> In this particular case we can try to see what happens when we use
> a variable limit in place of 2 Pi
>
> f = (1 + k*Sin[t]^2)^(1/2);
> q = Integrate[f, {t, 0, b}]
> If[k*Cos[2*b] != k + 2 && (Im[k*Sin[b]^2] != 0 ||
>     Re[k*Sin[b]^2] + 1 >= 0) &&
>    (Re[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] >= 0 ||
>     Im[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] != 0 ||
>     (Re[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] + 2 <= 0 &&
>      (Re[(((-Cos[2*b])*k + k + 2)*Csc[b]^2)/k] >= 2 ||
>       Im[(((-Cos[2*b])*k + k + 2)*Csc[b]^2)/k] != 0))),
>   EllipticE[b, -k], Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0, b},
>    Assumptions -> k*Cos[2*b] == k + 2 || (Im[k*Sin[b]^2] == 0 &&
>       Re[k*Sin[b]^2] + 1 < 0) ||
>      (Im[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] == 0 &&
>       Re[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] < 0 &&
>       (Re[((Cos[2*b]*k - k - 2)*Csc[b]^2)/k] + 2 > 0 ||
>        (Im[(((-Cos[2*b])*k + k + 2)*Csc[b]^2)/k] == 0 &&
>         Re[(((-Cos[2*b])*k + k + 2)*Csc[b]^2)/k] < 2)))]]
>
> This looks complicated, but
>
> Simplify[q /. b -> 2*Pi]
> 4*EllipticE[-k]
>
> which suggests that the result that you got under the restriction
> k >= -1 is actually valid without any restrictions on k. Empirical
> tests, with various values of k (real and complex) seem to confirm
> this.
>
> On the other hand, however:
>
>  Integrate[f, {t, 0, 2*Pi}]
>
> If[Re[k] > -1, 2*(EllipticE[-k] + Sqrt[k + 1]*EllipticE[k/(k + 1)]),
>   Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0, 2*Pi},
>    Assumptions -> Re[k] <= -1]]
>
> The problem is not only the we now have again a condition on Re[k]
> but also that this answer does not agree with the answer 4*EllipticE
> [-k] that we got earlier and, in particular,  with:
>
> Assuming[Element[k, Reals], Integrate[f, {t, 0, 2*Pi}]]
>
> If[k >= -1, 4*EllipticE[-k], Integrate[Sqrt[k*Sin[t]^2 + 1], {t, 0,
> 2*Pi}, Assumptions -> k < -1]]
>
> For example, taking k==5 (so that k>=-1) we get:
>
>  N[2*(EllipticE[-k] + Sqrt[k + 1]*EllipticE[k/(k + 1)]) /. k -> 5, 10]
>  11.32079299
>
> N[4 EllipticK[-k] /. k -> 5, 10]
> 3.822015708
>
> Andrzej Kozlowski
>

```

• Prev by Date: Re: [Mathematica 6] Integrate strange result
• Next by Date: Re: Suggestions for Maintaining "Object" State?
• Previous by thread: Re: [Mathematica 6] Integrate strange result
• Next by thread: Re: [Mathematica 6] Integrate strange result