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