Mathematica 9 is now available
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: [mg82173] Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Sun, 14 Oct 2007 06:05:14 -0400 (EDT)
  • References: <20071012203248.805$fe_-_@newsreader.com> <ACC57C97-2951-4DE6-829F-4E31B2BDFDF8@mimuw.edu.pl>

I still feel sufficiently provoked by your post to write one more  
reply, not to change what I wrote earlier in any fundamental way but  
just to make absolutely clear what the point of my first reply (to  
Bobby Treat) was, since it it appears to me that you have entirely  
missed it.


I will start with my first response to Bobby Treat.  Originally I did  
not intend to join this thread (since I consider this topic to be a  
horse that was beaten to death on this list a long time ago mostly by  
yourself). But in his (otherwise correct) reply, Bobby Treat  
expressed doubts about the correctness of the functions Simplify or D  
(yes, it was actually D rather than Integrate). Now, I know Simplify  
to be very reliable; there might have been perhaps one or two cases  
when it did not work correctly in the whole history of Mathematica.  
As for D, I would be really shocked to see it return a wrong answer,  
simply because it uses such a simple algorithm. I don't think it has  
ever done so. So that claim made interested in this thread and as I  
wrote, I immediately begin to suspect that the Newton-Leibniz rule  
was being used incorrectly as this is quite common not just in  
Mathematica but in any program that performs symbolic definite  
integration. However, while writing my reply I forgot that the it was  
D about which doubt had been expressed and mentioned Integrate.

Now, while Integrate is quite reliable it is not nearly as reliable  
as D, for reasons that I will try to explain below. But it is still  
very reliable in doing what it is supposed to do: finding "indefinite  
integrals". However, it is important to understand what is meant by  
an indefinite integral in symbolic algebra. The concept of an anti- 
derivative used in symbolic algebra is purely algebraic and does not  
depend on the idea of limit or on anything taught in undergraduate  
calculus. The whole process works over an arbitrary differential  
fields in the sense of Ritt. Differential fields are just fields with  
a D operator, that satisfies properties analogous to usual  
differentiation. If has a subfield of constants k such that D[k]=0,  
and one can consider algebraic and transcendental extensions in the  
usual sense. The problem of finding an "antiderivative" is the  
following: given an element f of a given differential field F, find  
an "allowable" differential extension G = F[x1,x2,...] and an element  
g in G such that D[g]=f or determine that no such "allowable"  
extension exists. If we take F to be just the field of rational  
functions over the rationals then there are several algorithms, which  
solve this problem and the extensions that are needed just usual  
algebraic field extensions of the rationals and transcendental  
(actually logarithmic) extensions of these. Of course when you think  
of these objects as the familiar objects from calculus they turn out  
to have branch cuts and singularities. The Risch algorithm of 1968  
solves the algebraic integration problem for a class of functions  
known as "elementary functions" (the existence of such a solution was  
proved much earlier by Liouville). It is exactly what this algorithm  
computes that is called an "indefinite integral" or an "anti- 
derivative" in the context of symbolic algebra systems.

The problem with the Risch algorithm is that it has two different  
aspects. The easy one is when we are able to confine ourselves to the  
class of purely transcendental extensions (which turn out to be just  
logarithmic and exponential). The hard one is when one needs to  
compute also algebraic extensions. As far as I know no symbolic  
algebra program implements completely the latter part. Instead some  
additionalo heuristics is used at some point or perhaps the  
integratal is returned unevaluated (?) This is why I would not  
consider Integrate to be nearly as reliable as D.

I do not know exactly what exactly is the procedure used by  
Mathematica's Integrate, but an analogous one used by one of its  
rivals is described in the book of Geddes, Czapor and Labahn which I  
mentioned in my earlier post. That is, given a function f whose anti- 
derivative is supposed to be found, the program first tries a number  
of heuristic methods, similar to those taught in calculus courses. It  
also makes use of some look-up tables. When all this fails, the  
program turns to the Risch algorithm. I am pretty sure Mathematica  
behaves in the same way.

Nowhere, at any stage of this process does any consideration of  
limits or  singularities or branch cuts come into play. These things  
have nothing to do with the purely algebraic process of indefinite  
integration. Integrate produces a "correct answer" unless there is a  
bug somewhere in its heuristics. Another possibility of course is  
that the Risch algorithm will take too long and the user aborts.

So in view of the above I can't make any sense of you claim that  
Integrate should be doubted or is not reliable unless what you  
actually meant was that using the anti-derivative returned by
Integrate and the  Newton-Leibniz formula is not reliable. This,
however, is just a  truism and in fact exactly what I wrote in my
post. Anyway, this  whole belongs to the realm of definite and not indefinite  
integration. Daniel Lichtblau has already posteda reply to Bobby the  
same thread (have you looked at it?)  where he included a link to a  
paper of his which explains in detail the difficulties faced by  
definite integration in this kind of setting. Definite integration  
still makes use of the Risch algorithm (and hence of derivatives with  
possible singular points) but it does not simple mindedly use the  
Newton-Leibniz formula. In fact, with numerical limits it correctly  
deals with the integral posted by the OP. It is only when the limits  
are symbolic that it fails to determine correctly the location of the  
branch cuts. Of course if it could use a purely real formula like  
yours in the case when definite integrals over an interval are  
involved, things would become much simpler i9n this particualr case,  
but this approach does not appear to be algorithmic enough to  
incorporate it into a general integration procedure. The only way to  
employ it would seem to be by means of a look-up table of definite  
integrals, which would deal with some cases but would certainly still  
leave (probably) infinitely many problematic ones.

Andrzej Kozlowski



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

>
> On 13 Oct 2007, at 09:32, David W. Cantrell wrote:
>
>> [Message also posted to: comp.soft-sys.math.mathematica]
>>
>> 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: Programming Euler's Method into Mathematica
  • Next by Date: format mixed integers & floats with text styling (see )
  • 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?