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: [mg82169] Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Sat, 13 Oct 2007 04:04:17 -0400 (EDT)
  • References: <20071012203248.805$fe_-_@newsreader.com> <ACC57C97-2951-4DE6-829F-4E31B2BDFDF8@mimuw.edu.pl>

Because of all this hurry I even managed to get the second  
fundamental theorem of calculus wrong. Of course it states that a  
continuous function on an interval has a C1-anti-derivative on that  
interval.

Andrzej Kozlowski


On 13 Oct 2007, at 10:23, Andrzej Kozlowski wrote:

>
> On 13 Oct 2007, at 09:32, David W. Cantrell wrote:
>
>>
>> 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.
>
> This has been discussed many time before on this list (with your  
> participation) but since you seem not to be satisfied, here it goes  
> again, I hope for the last time.
>
> The second fundamental theorem of calculus i.e. the Newton-Leibniz  
> rule (high school level mathematics) states that if we have a  
> function with a continuous derivative (a C1 function) on an  
> interval than it has a continuous anti-derivative on this  
> interval.   However if the original function is not complex  
> analytic in the complex plane then there will not be a complex  
> analytic anti-derivative. Instead, there will be "anti-derivatives"  
> with singularities, and in unlucky cases they may fall on the  
> interval of integration. This is exactly what happens 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.
>
> It depends on the meaning of "anti-derivative". In the meaning that  
> Mathematica and all other symbolic algebra programs use, a function  
> can be an an anti-derivative of another function, even if it is  
> singular or undefined  at certain points in the complex plane  
> provided it is an anti-derivative at almost all points (i.e. a  
> generic anti-derivative). It is of course completely natural for  
> Mathematica to use such anti-derivatives since it is precisely this  
> kind of anti-derivatives that are found by the Risch algorithm.
>
>>
>>> 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.
>
> Mathematica does not even attempt to find anti-derivatives on R in  
> your sense. Instead of bringing up time and again same points that  
> are covered in undergraduate calculus course I suggest learning a  
> little about symbolic algebraic integration. Here are two good  
> references:
>
> Geddes, Czapor, Labhan "Algorithms for Computer Algebra", Kluwer, 1992
> Davenport J. "On the Integration of Algebraic Functions", Springer- 
> Verlag 1981.
>
> If you read one of these you will find that the concept of "anti- 
> derivative" that your are using is not considered suitable for  
> computer algebra  (at least at its present stage of development  
> unless there have been some dramatic developments recently that I  
> am not aware of ). So, "you do not know how to obtain your answer  
> with Mathematica" because it can't be done. Symbolic algebra  
> programs use known algebraic algorithms. There is no point  
> demanding that they should do something unless a suitable general  
> algorithm is known.
> Of course, it is possible that you know an algorithmic (in the same  
> sense that the Risch algorithm is "algorithmic") way of computing  
> such anti-derivatives. If so, I am sure that Wolfram Research will  
> be very keen to implement it.
> But until then this is all there is to this entire issue.
>
> Andrzej Kozlowski
>
>
>
>>
>> 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: Can Integrate[expr,{x,a,b}] give an incorrect result?
  • Next by Date: Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
  • Previous by thread: 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?