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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg81912] Re: [mg81855] Re: [mg81827] Can Integrate[expr, {x, a, b}] give an
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Sat, 6 Oct 2007 04:44:29 -0400 (EDT)
  • References: <8207246.1191505119218.JavaMail.root@m35> <200710050846.EAA00304@smc.vnet.net>

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 > 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
>>

To summarize: Mathematica has computed the antiderivative, evaluated at 
the endpoints, subtracted, and obtained zero. That zero is the correct 
difference of antiderivatives at endpoints of integration. Moreover, as 
Bobby Treat points out, there is a branch cut crossed, and Integrate 
fails to account for it by splitting the integration path. Neither D not 
Simplify has made a mistake above, at least not that I could discern.

GenerateConditions default behavior, for univariate integration, is to 
try to generate them. So setting it explicitly to True is unlikely to 
elicit better behavior for this sort of problem. It is already trying, 
and not succeeding (or, more bluntly, it's failing).

As for the missed branch cut crossing, this seems to be a weakness in 
handling parametrized solution sets. For alpha>0 there are (countably) 
infinitely many values for which this will happen. The code utilized to 
detect them will at best find finitely many (possibly only those given 
by Solve in arctrig solutions, i.e. between -Pi and Pi or 0 and 2*Pi). A 
way to force a better result, for this problem, is to restrict alpha to, 
say, 0<alpha<2*Pi.

Apologies if you spent more time on this than you wanted to. If it is 
any consolation, so did I. That is to say, I'm paid to look at this 
stuff, but it still sometimes hurts to do so.


Daniel Lichtblau
Wolfram Research





  • Prev by Date: Re: Re: converting function from RasterArray to Raster
  • Next by Date: Re: Help notebook style in Mathematica 6.0
  • Previous by thread: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
  • Next by thread: Re: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?