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

MathGroup Archive 2007

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

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg81855] Re: [mg81827] Can Integrate[expr,{x,a,b}] give an incorrect result?
  • From: DrMajorBob <drmajorbob at bigfoot.com>
  • Date: Fri, 5 Oct 2007 04:46:36 -0400 (EDT)
  • References: <8207246.1191505119218.JavaMail.root@m35>
  • Reply-to: drmajorbob at bigfoot.com

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



-- 

DrMajorBob at bigfoot.com


  • Prev by Date: Re: Re: Problems with PlotStyle in ListPointPlot3D
  • Next by Date: Re: Manipulate+Plot showing no output
  • Previous by thread: Result Duplicated! Re: Can Integrate[expr,{x,a,b}] give
  • Next by thread: Re: Re: Can Integrate[expr, {x, a, b}] give an