Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
- To: mathgroup at smc.vnet.net
- Subject: [mg82158] Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
- From: "David W.Cantrell" <DWCantrell at sigmaxi.net>
- Date: Sat, 13 Oct 2007 03:58:36 -0400 (EDT)
- References: <8207246.1191505119218.JavaMail.root@m35> <200710050846.EAA00304@smc.vnet.net> <fe7h33$od3$1@smc.vnet.net>
Sorry for not having responded sooner; I've been away from computer access
for a week.
Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote:
> > 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).
Of course, the Leibniz rule cannot be used unless one _has_ an
antiderivative valid on the required interval. Mathematica does not supply
such an antiderivative in this case.
> 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?
It _must_ give the right answer if indef is actually an antiderivative on
the interval 0 to 2Pi, but it isn't.
> The function indef is clearly discontinuous in the interval 0 to 2Pi
True, but that discontinuity is merely an artifact of Mathematica's
workings. The integrand int is clearly continuous on the interval 0 to 2Pi
and, indeed, on all of R, and so there is an antiderivative on R:
886400 Sqrt[377]/4340193273 (t - 2 ArcTan[44 Sin[t]/(44 Cos[t] - 3 Sqrt[377] - 73)]) + 8800 Sin[t] (1078 Cos[t] - 8009)/(3837483 (44 Cos[t] - 73)^2))
but I don't know how to get Mathematica to produce it. Of course, using my
antiderivative above with Leibniz rule, we get the correct result for the
definite integral.
> 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.
There is very good reason to suspect that Integrate, when asked to evaluate
an indefinite integral, might give a result which is not an antiderivative
on some desired interval. This happens fairly often, the subject of this
thread being a case in point.
Another result which might help Craig, but which I don't know how to get
using Mathematica, is that
Assuming[r > 0 && rho > 0 && Element[zeta, Reals] && Element[z, Reals] &&
alpha > 0 && Not[rho == r && zeta == z], Integrate[(r*(r - rho*Cos[alpha -
t]))/ (3*(r^2 - 2*rho*Cos[alpha - t]*r + rho^2 + (z - zeta)^2)^3), {t, 0, 2 Pi}]]
should be
(2*Pi*r^2*(-2*rho^4 + rho^2*(r^2 - (zeta - z)^2) + (r^2 + (zeta - z)^2)^2))/
(3*(((rho - r)^2 + (zeta - z)^2)*((rho + r)^2 + (zeta - z)^2))^(5/2))
Note that the last assumption mentioned above was added by me. Without that
last assumption, there would have been trouble if we had had both zeta == z
and rho == r. But in that special case, the desired definite integral is
easily handled: it's just 0.
David W. Cantrell
> 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
> >
- 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?