Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
- To: mathgroup at smc.vnet.net
- Subject: [mg82168] 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:03:46 -0400 (EDT)
- References: <20071012203248.805$fe_-_@newsreader.com> <ACC57C97-2951-4DE6-829F-4E31B2BDFDF8@mimuw.edu.pl>
Just to avoid any more disputes ove relementary calculus, a minor
correction to the following (which was written in a hurry):
> 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.
Actually, the Newton-Leibniz formula is "The first fundamental
theorem of Calculus" and the existence of a continuos anti-derivative
of a C1-fucntion on an interval follows from "The second fundamental
theorem of Calculus".
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
>>>>
>