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