Re: Can Integrate[expr,{x,a,b}] give an incorrect result?
- To: mathgroup at smc.vnet.net
- Subject: [mg82165] 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:02:12 -0400 (EDT)
- References: <20071012203248.805$fe_-_@newsreader.com>
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
>>>