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
- References:
- Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
- From: DrMajorBob <drmajorbob@bigfoot.com>
- Re: Can Integrate[expr,{x,a,b}] give an incorrect result?