Re: Definite Integration in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg74498] Re: Definite Integration in Mathematica*From*: "dimitris" <dimmechan at yahoo.com>*Date*: Fri, 23 Mar 2007 19:07:08 -0500 (EST)*References*: <etqo3f$10i$1@smc.vnet.net><ett73k$h3g$1@smc.vnet.net>

You miss one important thing! (x^2 + 2*x + 4)/(x^4 - 7*x^2 + 2*x + 17) (the function in discussion) is indeed continuous in the whole real axis. (We don't talk just for Reimmanian integration in the real axis as most of calculus books deal with! It is rather an issue of Complex Analysis.) But what about the complex plane? If it was analytic everywhere in the complex plane then it would have an analytic anti- derivative everywhere, but here is not the case. The function has four poles in the complex plane! Also convince yourself that the antiderivative returned by Mathematica is indeed right. (there is of course, in general, a whole family of results which have the property that their derivative is the integrand). In[237]:= Derivative[-1][(#1^2 + 2*#1 + 4)/(#1^4 - 7*#1^2 + 2*#1 + 17) & ][x] Simplify[%] D[%, x] Simplify[%] Out[237]= (1/2)*ArcTan[(-1 - x)/(-4 + x^2)] - (1/2)*ArcTan[(1 + x)/(-4 + x^2)] Out[238]= ArcTan[(1 + x)/(4 - x^2)] Out[239]= ((2*x*(1 + x))/(4 - x^2)^2 + 1/(4 - x^2))/(1 + (1 + x)^2/(4 - x^2)^2) Out[240]= (4 + 2*x + x^2)/(17 + 2*x - 7*x^2 + x^4) In another symbolic system (again the policy of MathGroup does not allow to mention its name) I get indeed a continuous antiderivative in the whole real axis int((x^2 + 2*x + 4)/(x^4 - 7*x^2 + 2*x + 17),x); arctan(5/3-x-1/3*x^2+1/3*x^3)+arctan(x-1) plot(%,x=-10..10); but this does not say anything for the integration algorithm of Mathematica! As the function is not analytic it will only have an analytic anti- derivative in any simply connected region in which it is analytic. So any candidate for an antidrivative of this function, will have to have singularities somewhere in the region. Before proceeding a check is done In[246]:= D[other, x] Simplify[%] Out[246]= 1/(1 + (1 - x)^2) + (-1 - (2*x)/3 + x^2)/(1 + (5/3 - x - x^2/3 + x^3/3)^2) Out[247]= (4 + 2*x + x^2)/(17 + 2*x - 7*x^2 + x^4) Let compare now the two candidate antiderivatives In[250]:= math = ArcTan[(1 + x)/(4 - x^2)]; (*mathematica's antiderivative*) other = ArcTan[5/3 - x - (1/3)*x^2 + (1/3)*x^3] + ArcTan[x - 1]; (*other CAS antiderivative*) See the following graphs In[254]:= Show[GraphicsArray[Block[{$DisplayFunction = Identity}, (ContourPlot[Evaluate[#1[math /. x -> x + I*y]], {x, -4, 4}, {y, -3, 3}, ContourShading -> False, PlotPoints -> 50, Contours -> 50] & ) /@ {Re, Im}], Frame -> True], ImageSize -> 600, PlotLabel -> "math antiderivative"]; Show[GraphicsArray[Block[{$DisplayFunction = Identity}, (ContourPlot[Evaluate[#1[other /. x -> x + I*y]], {x, -4, 4}, {y, -3, 3}, ContourShading -> False, PlotPoints -> 50, Contours -> 50] & ) /@ {Re, Im}], Frame -> True], ImageSize -> 600, PlotLabel -> "other CAS antiderivative"]; Show[GraphicsArray[Block[{$DisplayFunction = Identity}, (Plot3D[Evaluate[#1[math /. x -> x + I*y]], {x, -4, 4}, {y, -3, 3}, PlotPoints -> 50] & ) /@ {Re, Im}], Frame -> True], ImageSize -> 800, PlotLabel -> "math antiderivative"]; Show[GraphicsArray[Block[{$DisplayFunction = Identity}, (Plot3D[Evaluate[#1[other /. x -> x + I*y]], {x, -4, 4}, {y, -3, 3}, PlotPoints -> 50] & ) /@ {Re, Im}], Frame -> True], ImageSize -> 800, PlotLabel -> "other CAS antiderivative"]; I think you are convinced now that both antiderivatives have discontinuities in the complex plane! It's just happens the returned antiderivative from the other CAS to be continuous in the real axis; it's a matter of choice of the region of analyticity! Note also that In[260]:= math - other Plot[Chop[%], {x, -3, 3}, Axes -> False, Frame -> True] I=2Ee. in the real axis the two antiderivatives differ by a piecewise (complex as we saw in the next plot) constant function and inside (-2,2) are equally! In[264]:= Show[GraphicsArray[Block[{$DisplayFunction = Identity}, (Plot3D[Chop[Evaluate[#1[math - other /. x -> x + I*y]]], {x, -4, 4}, {y, -3, 3}, PlotPoints -> 50] & ) /@ {Re, Im}], Frame -> True], ImageSize -> 800]; Obviously In[271]:= D[math - other, x] Simplify[%] Out[271]= -(1/(1 + (1 - x)^2)) + ((2*x*(1 + x))/(4 - x^2)^2 + 1/(4 - x^2))/(1 + (1 + x)^2/(4 - x^2)^2) - (-1 - (2*x)/3 + x^2)/(1 + (5/3 - x - x^2/3 + x^3/3)^2) Out[272]= 0 So everything is normal regarding the application of the Newton Leibniz formula as the following demonstrates: In[279]:= FullSimplify[(math /. x -> 4) - Limit[math, x -> 2, Direction -> -1] + Limit[math, x -> 2, Direction -> 1] - math /. x -> 0] N[%] FullSimplify[(other /. x -> 4) - (other /. x -> 0)] N[%] Out[279]= Pi - (1/2)*ArcTan[2752/825] Out[280]= 2=2E5018228707631676 Out[281]= (1/4)*(3*Pi + 2*ArcTan[825/2752]) Out[282]= 2=2E5018228707631676 Then the question is: Can we use Mathematica to get an antiderivative continuous in the real axis? Certainly we can for this example. But usually the process of builting an antiderivative continuous in a given integral is a very tough procedure. But, note that we can't get rid of the brach cuts! (At last I think I understand Adamchik's article; Thanks a lot Daniel Huber and Andrzej Kozlowski!) So for the our integrand (copy/paste the following as a WHOLE cell) In[1471]:= Print["the integrand"] f[x_] = (4 + 2*x + x^2)/(17 + 2*x - 7*x^2 + x^4) Print["\[Cross]partial fraction decomposition procedure\[Cross]"] Print["the denominator of f[x] as a product of terms"] Times @@ Apply[#1[[1]] - #1[[2]] & , Solve[Denominator[f[x]] == 0, x], 1] Print["f[x] as a sum of fractions"] Map[FullSimplify, Apart[(4 + 2*x + x^2)/%%], 1] Print["Integrate each term of above sum"] F = (Integrate[#1, x] & ) /@ %% Print["plot of the antiderivative"] Plot[%%, {x, -10, 10}, Axes -> False, Frame -> True]; Print["verification"] FullSimplify[D[F, x] == f[x]] Print["the definite integral in the interval [0,4]"] (F /. x -> 4) - (F /. x -> 0) Print["numerical value"] %%//N//Chop Do you prefer the extend antiderivative over the compact one obtained directly by Mathematica only because is real in the real axis? As regards myself, no! Of course the question of simplifications arises now but this is not then subject of the current thread! Regards Dimitris /=C7 Michael Weyrauch =DD=E3=F1=E1=F8=E5: > Hello, > > another nice example, where the result for the integral given by Mathem= atica just > cannot be right. The indefinite integral of a continuuos function cannot = have a jump. > That -- to my opinion -- is mathematical nonsens. We all learned that int= egration "smoothens", > i.e. if the integrand is somewhat "ugly" the integral is less "ugly". (Ne= ver mind my English!). > > So the result should actually be presented as > > F[x_]=ArcTan[(1 + x)/(4 - x^2)]*UnitStep[2 - x] + > (Pi + ArcTan[(1 + x)/(4 - x^2)])*UnitStep[-2 + x] > > This is a perfectly nice function without jumps and it is the antideriva= tive of your integrand, > and the fundamental theorem of calculus works with it... > > So, despite my love for Mathematica, here it fools me.... > > Regards Michael > > > > > "dimitris" <dimmechan at yahoo.com> schrieb im Newsbeitrag news:etqo3f$10i$1= @smc.vnet.net... > > Hello to all of you! > > > > Firstly, I apologize for the lengthy post! > > Secondly, this post has a close connection with a recent (and well > > active!) > > thread titled "Integrate" and one old post of mine which was based on > > a older > > post of David Cantrell. Since there was no response and I do consider > > the > > subject very fundamental I would like any kind of insight. > > > > In the section about Proper Integrals in his article, Adamchik > > mentions that the Newton-Leibniz formula (i.e. the Fundamental > > Theorem of Integral Calculus: Integrate[f[x],{x,a,b}]=F[b]-F[a], > > F[x]: an antiderivative), does not hold any longer if the > > antiderivative F(x) has singularities in the integration interval > > (a,b). > > > > To demonstrate this, he considers the integral of the function: > > > > f[x_] = (x^2 + 2*x + 4)/(x^4 - 7*x^2 + 2*x + 17); > > > > over the interval (0,4). > > > > Plot[f[x], {x, 0, 4}]; > > (*plot to be displayed*) > > > > The integrand posseses no singularities on the interval (0,4). > > > > Here is the corresponding indefinite integral > > > > F[x_] = Simplify[Integrate[f[x], x]] > > ArcTan[(1 + x)/(4 - x^2)] > > > > Substituting limits of integration into F[x] yields an incorrect > > result > > > > Limit[F[x], x -> 4, Direction -> 1] - Limit[F[x], x -> 0, Direction - > >>1] > > N[%] > > NIntegrate[f[x], {x, 0, 4}] > > > > -ArcTan[1/4] - ArcTan[5/12] > > -0.6397697828266257 > > 2.501822870767894 > > > > This is because the antiderivative has a jump discontinuity at x=2 > > (also at x = -2), so that the Fundamental theorem cannot be used. > > > > Indeed > > > > Limit[F[x], x -> 2, Direction -> #1]&/@{-1, 1} > > Show@Block[{$DisplayFunction=Identity}, > > Plot[F[x],{x,#[[1]],#[[2]]}]&/@Partition[Range[0,4,2],2,1]]; > > {-(Pi/2), Pi/2} > > (*plot to be displayed*) > > > > The right way of applying the Fundamental theorem is the following > > > > (Limit[F[x], x -> 4, Direction -> 1] - Limit[F[x], x -> 2, Direction - > >> -1]) + > > (Limit[F[x], x -> 2, Direction -> 1] - Limit[F[x], x -> 0, Direction > > -> -1]) > > N[%] > > Pi - ArcTan[1/4] - ArcTan[5/12] > > 2.501822870763167 > > > > Integrate works in precisely this way > > > > Integrate[f[x], {x, 0, 4}] > > N[%] > > > > Pi - ArcTan[1/4] - ArcTan[5/12] > > 2.501822870763167 > > > > A little later, he (i.e. Adamchik) says "The origin of > > discontinuities > > along the path of integration is not in the method of indefinite > > integration but rather in the integrand." > > > > Adamchik mentions next that the four zeros of the integrand's > > denominator > > are two complex-conjugate pairs having real parts +/- 1.95334. It > > then > > seems that he is saying that, connecting these conjugate pairs by > > vertical line segments in the complex plane, we get two branch > > cuts... > > > > BUT didn't the relevant branch cuts for his int cross the real axis > > at x = +/- 2, rather than at x = +/- 1.95334? > > > > (NOTE: The difference between 1.95334 and 2 is not due to numerical > > error). > > > > Exactly what's going on here? > > > > Show[GraphicsArray[Block[{$DisplayFunction = Identity}, > > (ContourPlot[#1[F[x + I*y]], {x, -4, 4}, {y, -4, 4}, Contours -> 50, > > PlotPoints -> 50, ContourShading -> False, Epilog -> {Blue, > > AbsoluteThickness[2], Line[{{0, 0}, {4, 0}}]}, > > PlotLabel -> #1[HoldForm[F[x]]]] & ) /@ {Re, Im}]], ImageSize > > - > 500]; > > (*contour plots to be displayed*) > > > > > > Consider next the following function > > > > f[x_] = 1/(5 + Cos[x]); > > > > Then > > > > Integrate[f[x], {x, 0, 4*Pi}] > > N[%] > > NIntegrate[f[x], {x, 0, 4*Pi}] > > > > Sqrt[2/3]*Pi > > 2.565099660323728 > > 2.5650996603270704 > > > > F[x_] = Integrate[f[x], x] > > ArcTan[Sqrt[2/3]*Tan[x/2]]/Sqrt[6] > > > > D[F[x], x]==f[x]//Simplify > > True > > > > Plot[f[x], {x, 0, 4*Pi}, Ticks -> {Range[0, 4*Pi, Pi/2], Automatic}] > > Plot[F[x], {x, 0, 4*Pi}, Ticks -> {Range[0, 4*Pi, Pi/2], Automatic}] > > > > The antiderivative has jump discontinuities at Pi and 3Pi inside the > > integration range > > > > Table[(Limit[F[x], x -> n*(Pi/2), Direction -> #1] & ) /@ {-1, 1}, {n, > > 0, 4}] > > {{0, 0}, {ArcTan[Sqrt[2/3]]/Sqrt[6], ArcTan[Sqrt[2/3]]/Sqrt[6]}, {-(Pi/ > > (2*Sqrt[6])), Pi/(2*Sqrt[6])}, > > {-(ArcTan[Sqrt[2/3]]/Sqrt[6]), -(ArcTan[Sqrt[2/3]]/Sqrt[6])}, {0, > > 0}} > > > > Reduce[5 + Cos[x] == 0 && 0 <= Re[x] <= 4*Pi, x] > > {ToRules[%]} /. (x_ -> b_) :> x -> ComplexExpand[b] > > x /. %; > > ({Re[#1], Im[#1]} & ) /@ %; > > poi = Point /@ %; > > > > x == 2*Pi - ArcCos[-5] || x == 4*Pi - ArcCos[-5] || x == Ar= cCos[-5] || > > x == 2*Pi + ArcCos[-5] > > {{x -> Pi - I*Log[5 - 2*Sqrt[6]]}, {x -> 3*Pi - I*Log[5 - 2*Sqrt[6]]}, > > {x -> Pi + I*Log[5 - 2*Sqrt[6]]}, > > {x -> 3*Pi + I*Log[5 - 2*Sqrt[6]]}} > > > > Of course F[4Pi]-F[0]=0 incorrectly. > > > > The reason for the discrepancy in the above result is not because of > > any problem with the fundamental theorem of calculus, of course; it is > > caused by the multivalued nature of the indefinite integral arctan. > > > > > > Show[GraphicsArray[Block[{$DisplayFunction = Identity}, > > (ContourPlot[#1[F[x + I*y]], {x, 0, 4*Pi}, {y, -4, 4}, Contours -> > > 50, PlotPoints -> 50, ContourShading -> False, > > FrameTicks -> {Range[0, 4*Pi, Pi], Automatic, None, None}, > > Epilog -> {{PointSize[0.03], Red, poi}, > > {Blue, Line[{{0, 0}, {4*Pi, 0}}]}}] & ) /@ {Re, Im}]], > > ImageSize -> 500] > > > > So, in this example the discontinuities are indeed from the branch > > cuts that start and end from the simple poles of the integrand which > > is in agreement with V.A. paper! > > > > > > I think I am not aware of something fundamental! > > Can someone point out what I miss? > > > > > > Regards > > Dimitris > > > >