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

```

• Prev by Date: Re: Re: [Mathematica 6] Format->Magnification does not work on notebook
• Next by Date: Re: How to do count for sub list?????
• Previous by thread: Re: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
• Next by thread: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?