Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

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

Search the Archive

Re: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg81883] Re: [mg81855] Re: [mg81827] Can Integrate[expr,{x,a,b}] give an incorrect result?
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Sat, 6 Oct 2007 04:29:33 -0400 (EDT)
  • References: <8207246.1191505119218.JavaMail.root@m35> <200710050846.EAA00304@smc.vnet.net>

> So, from that, D and/or Simplify must be wrong.

Must they? Even without any investigation I would say that it is much  
more likely that this integral should not be evaluated by means of  
the Leibniz rule (substituting the limits into the anti-derivative  
and subtracting). In fact, it is easy to see what happens in your  
special case:

integrand = (r*(r - rho*Cos[alpha - t]))/
        (3*(r^2 - 2*rho*Cos[alpha - t]*r + rho^2 +
                (z - zeta)^2)^3);
  int = integrand /. {r -> 1, rho -> 11/10,
        zeta -> 0, alpha -> 0, z -> 12/10}
  (1 - (11*Cos[t])/10)/(3*(73/20 - (11*Cos[t])/5)^3)

First, let us compute the integral using Integrate:

Integrate[int, {t, 0, 2*Pi}]
(1772800*Pi)/(11512449*Sqrt[377])

N[%]
0.0249156

Now, let's do it by numerical integration:

  NIntegrate[int, {t, 0, 2*Pi}]
  0.024915592618816908

so far, so good. Now let's compute the indefinte integral:

indef = Integrate[int, t];

Using the Liebniz rule gives clearly the wrong answer:

indef = Integrate[int, t];

but why should it give the right answer? The function indef is  
clearly discontinuous in the interval 0 to 2Pi

Plot[indef, {t, 0, 2 Pi}]

so the Leibniz rule does not apply. There is no reason at all to  
suspect the very reliable function Simplify and almost as reliabel  
Interate.

Andrzej Kozlowski





Andrzej Kozlowski



On 5 Oct 2007, at 17:46, DrMajorBob wrote:

> Yes, it's an incorrect result.
>
> To paraphrase your post:
>
> integrand = (r*(r -
>         rho*Cos[alpha - t]))/(3*(r^2 + rho^2 + (z - zeta)^2 -
>          2*r*rho*Cos[alpha - t])^3);
> assumptions = {r > 0, rho > 0, Element[zeta, Reals],
>     Element[z, Reals], alpha >In[15]:= 0};
>
> Integrate[integrand, {t, 0, 2 Pi}, Assumptions -> assumptions]
>
> 0
>
> Plot[integrand /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 0,
>     z -> 1.2}, {t, 0, 2 Pi}, PlotRange -> All]
>
> "Visual integration" is unconvincing, but you're right, the  
> integral isn't
> zero:
>
> NIntegrate[
>   integrand /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 0,
>     z -> 1.2}, {t, 0, 2 Pi}]
>
> 0.0249156
>
> Your example violates "assumptions" (alpha -> 0 vs. alpha > 0), but  
> that's
> not the problem:
>
> NIntegrate[
>   integrand /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 1,
>     z -> 1.2}, {t, 0, 2 Pi}]
>
> 0.0249156
>
> So, as often happens, you've found a case where Mathematica got a  
> general
> solution that isn't correct for some values of the parameters. The
> quadratic formula is another example:
>
> Solve[a x^2 + b x + c == 0, x]
>
> {{x -> (-b - Sqrt[b^2 - 4 a c])/(
>     2 a)}, {x -> (-b + Sqrt[b^2 - 4 a c])/(2 a)}}
>
> Not true when a == 0, of course.
>
> But let's look further:
>
> indefinite = Integrate[integrand, t, Assumptions -> assumptions];
> definite = Subtract @@ (indefinite /. {{t -> 2 Pi}, {t -> 0}});
> Simplify[definite, assumptions]
>
> 0
>
> Is that a Simplify error, or an Integrate error?
>
> definite /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 1, z -> 1.2}
>
> -2.05998*10^-18 + 0. \[ImaginaryI]
>
> Ah. Simplify may be doing its job just fine (hard to be certain),  
> but it
> appears the indefinite integral is incorrect. But look:
>
> D[indefinite, t] == integrand // Simplify
>
> True
>
> and
>
> D[indefinite, t] == integrand // Simplify[#, assumptions] &
>
> True
>
> and also
>
> definite // Together // Numerator // Simplify
>
> 0
>
> So, from that, D and/or Simplify must be wrong.
>
> If this is correct:
>
> numdefinite // Together // Numerator // TrigFactor
>
> -4 r^2 (r^4 + r^2 rho^2 - 2 rho^4 + 2 r^2 z^2 - rho^2 z^2 + z^4 -
>     4 r^2 z zeta + 2 rho^2 z zeta - 4 z^3 zeta + 2 r^2 zeta^2 -
>     rho^2 zeta^2 + 6 z^2 zeta^2 - 4 z zeta^3 +
>     zeta^4) (ArcTanh[((r^2 + 2 r rho + rho^2 + (z - zeta)^2) Tan 
> [alpha/
>        2])/Sqrt[-r^4 +
>       2 r^2 (rho^2 - (z - zeta)^2) - (rho^2 + (z - zeta)^2)^2]] -
>     ArcTanh[((r^2 + 2 r rho + rho^2 + (z - zeta)^2) Tan[
>        1/2 (alpha - 2 \[Pi])])/
>      Sqrt[-r^4 +
>       2 r^2 (rho^2 - (z - zeta)^2) - (rho^2 + (z - zeta)^2)^2]])  
> (r^2 +
>      rho^2 + z^2 - 2 z zeta + zeta^2 - 2 r rho Cos[alpha])^2
>
> The next-to-last factor (the only one that depends on t) is
>
> num[[-2]]
>
> ArcTanh[((r^2 + 2 r rho + rho^2 + (z - zeta)^2) Tan[alpha/2])/
>    Sqrt[-r^4 +
>     2 r^2 (rho^2 - (z - zeta)^2) - (rho^2 + (z - zeta)^2)^2]] -
>   ArcTanh[((r^2 + 2 r rho + rho^2 + (z - zeta)^2) Tan[
>      1/2 (alpha - 2 \[Pi])])/
>    Sqrt[-r^4 + 2 r^2 (rho^2 - (z - zeta)^2) - (rho^2 + (z - zeta)^2) 
> ^2]
>    ]
>
> so it all comes down (seemingly) to
>
> Tan[alpha/2] == Tan[1/2 (alpha - 2 \[Pi])] // Simplify
>
> True
>
> or
>
> Tan[any] == Tan[any - Pi]
>
> True
>
> BUT... what's likely happening is that Simplify doesn't properly  
> take into
> account the branch-cut discontinuities of ArcTanh. In that case,  
> there's a
> significant set of parameters for which the integral IS zero, such as:
>
> Reduce[Flatten@{num == 0, assumptions}]
>
> ors = (zeta \[Element]
>      Reals && ((z < zeta &&
>         rho > Sqrt[z^2 - 2 z zeta + zeta^2]/Sqrt[2] && alpha > 0 &&
>         r == Root[-2 rho^4 - rho^2 z^2 + z^4 + 2 rho^2 z zeta -
>             4 z^3 zeta - rho^2 zeta^2 + 6 z^2 zeta^2 - 4 z zeta^3 +
>             zeta^4 + (rho^2 + 2 z^2 - 4 z zeta +
>                2 zeta^2) #1^2 + #1^4 &, 2]) || (z == zeta && rho >  
> 0 &&
>          alpha > 0 &&
>         r == Root[-2 rho^4 - rho^2 z^2 + z^4 + 2 rho^2 z zeta -
>             4 z^3 zeta - rho^2 zeta^2 + 6 z^2 zeta^2 - 4 z zeta^3 +
>             zeta^4 + (rho^2 + 2 z^2 - 4 z zeta +
>                2 zeta^2) #1^2 + #1^4 &, 2]) || (z > zeta &&
>         rho > Sqrt[z^2 - 2 z zeta + zeta^2]/Sqrt[2] && alpha > 0 &&
>         r == Root[-2 rho^4 - rho^2 z^2 + z^4 + 2 rho^2 z zeta -
>             4 z^3 zeta - rho^2 zeta^2 + 6 z^2 zeta^2 - 4 z zeta^3 +
>             zeta^4 + (rho^2 + 2 z^2 - 4 z zeta +
>                2 zeta^2) #1^2 + #1^4 &, 2]))) || ((z |
>        zeta) \[Element] Reals && alpha > 0 && r > 0 &&
>     rho > 0) || (C[1] \[Element] Integers && z \[Element] Reals &&
>     C[1] >= 1 && alpha == 2 \[Pi] C[1] && r > 0 && rho == r &&
>     zeta == z)
>
> But look:
>
> % /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 1, z -> 1.2}
>
> True
>
> Hold on, though... there are a lot of Ors in there. The least  
> restrictive
> seems to be
>
> ors[[1, 2, 1]]
>
> z < zeta && rho > Sqrt[z^2 - 2 z zeta + zeta^2]/Sqrt[2] && alpha >  
> 0 &&
>    r == Root[-2 rho^4 - rho^2 z^2 + z^4 + 2 rho^2 z zeta - 4 z^3  
> zeta -
>        rho^2 zeta^2 + 6 z^2 zeta^2 - 4 z zeta^3 +
>       zeta^4 + (rho^2 + 2 z^2 - 4 z zeta + 2 zeta^2) #1^2 + #1^4 &, 2]
>
> But that's a pretty strict condition (on r, especially), and
>
> % /. {r -> 1, rho -> 1.1, zeta -> 0, alpha -> 1, z -> 1.2}
>
> False
>
> SO. Integrate found a set of conditions for which the definite  
> integral
>
> was zero; it just didn't tell you what they were! This is supposed  
> to be
> the fix for that, I think:
>
> Integrate[integrand, {t, 0, 2 Pi}, Assumptions -> assumptions,
>   GenerateConditions -> True]
>
> 0
>
> BUT, as you see, that didn't help one bit!
>
> Bobby
>
> On Thu, 04 Oct 2007 03:27:21 -0500, W. Craig Carter <ccarter at mit.edu>
> wrote:
>
>>
>> I believe I am getting an incorrect result from a definite
>> integration:
>>
>> InputForm[integrand] is
>> (R*(R - rho*Cos[alpha - t]))/(3*(R^2 + rho^2 + (z - zeta)^2 -
>> 2*R*rho*Cos[alpha - t])^3)
>>
>> InputForm[assumptions] is
>> {R > 0, L > 0, rho > 0, Element[zeta, Reals], Element[z, Reals],  
>> alpha>
>> 0}
>>
>> Integrate[integrand,{t,0,2Pi},Assumptions->assumptions]
>> returns 0
>>
>> But compare this to:
>> (visually integrate...)
>> Plot[integrand/.{R -> 1, rho -> 1.1, zeta -> 0, alpha -> 0, z ->
>> 1.2},{t,0,2 Pi},PlotRange->All]
>>
>> (numerically integrate...)
>> Plot[NIntegrate[integrand/.{{R -> 1, rho -> 1.1, zeta -> 0, alpha - 
>> > 0,
>> z -> 1.2},{t,0,tau}],{tau,0,2Pi}, PlotRange->All]
>>
>>
>> Something isn't adding up??
>>
>> Thanks, WCC
>>
>>
>
>
>
> -- 
>
> DrMajorBob at bigfoot.com
>



  • Prev by Date: Re: Help notebook style in Mathematica 6.0
  • Next by Date: Re: What is the purpose of the Defer Command?
  • Previous by thread: Re: Re: Can Integrate[expr, {x, a, b}] give an
  • Next by thread: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?