Re: Integrate (a difficult integral)

*To*: mathgroup at smc.vnet.net*Subject*: [mg74515] Re: Integrate (a difficult integral)*From*: "dimitris" <dimmechan at yahoo.com>*Date*: Fri, 23 Mar 2007 19:16:28 -0500 (EST)*References*: <esls78$q0v$1@smc.vnet.net><et5o76$jor$1@smc.vnet.net>

As I said in another working with CAS don't mean that we... stop thinking with our mind. To obtain the best results (which are possible of course!)-not only in integration but in general doing serious Mathematics with Mathematica- you need an a combination of art, experience, reading, thinking (and guessing sometimes!). Let demonstrate one way of working for one integral that appeared in a recent thread in another forum. (The question was by Vladimir Bondarenko and I do consider that the integral should be appeared also in MathGroup.) Note that I don't say that it is the best one for the present situation but rather than it is a way that helped Mathematica to overcome its shothcomings in the present situation and return results in reasonable time (or someone would prefer to say: "Oh man, don't say many things; it's...just what you think! No more, no less!") Anyway... $Version "5.2 for Microsoft Windows (June 20, 2005)" Consider the function f[z_] := 1/Sqrt[z^3 - I*z^2 + z - I] The following integral returns unevaluated Integrate[F[z], {z, -Infinity, Infinity}] It is exactly this integral that we are looking for! The indefinite integral is F[z_] = Integrate[1/Sqrt[z^3 - I*z^2 + z - I], z] -(((1/4 + I/4)*(-I + z)*Sqrt[I + z]*(2*ArcTan[Sqrt[I + z]/(-2 + Sqrt[I + z])] - 2*ArcTan[Sqrt[I + z]/(2 + Sqrt[I + z])] + I*(Log[(2 + I) + z - 2*Sqrt[I + z]] - Log[(2 + I) + z + 2*Sqrt[I + z]])))/Sqrt[(-I + z)^2*(I + z)]) Of course "blinded" application of the Nweton-Leibniz formula leads to incorrect result Limit[F[z], z -> Infinity] - Limit[F[z], z -> -Infinity] {N[%], NIntegrate[f[x], {x, -Infinity, Infinity}]} (-(1/2) - I/2)*Pi {-1.5707963267948966 - 1.5707963267948966*I, 3.1415926536059606 + 3=2E14159265360596*I} Here are some plots (Plot[#1[f[z]], {z, -10, 10}] & ) /@ {Re, Im}; (Plot[#1[F[z]], {z, -20, 20}, Ticks -> {Range[-20, 20, 4], Automatic}, Axes -> False, Frame -> True] & ) /@ {Re, Im}; (ContourPlot[#1[F[x + I*y]], {x, -5, 5}, {y, -5, 5}, PlotPoints -> 100, Contours -> 100, ContourShading -> False] & ) /@ {Re, Im}; (ContourPlot[#1[f[x + I*y]], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 100, Contours -> 100, ContourShading -> False] & ) /@ {Re, Im} As we see both the real and imaginary part has jump discontinuities in the real axis. The following plots are very hepful (Plot[#1[ArcTan[Sqrt[I + z]/(-2 + Sqrt[I + z])]], {z, -10, 10}] & ) /@ {Re, Im}; (Plot[#1[ArcTan[Sqrt[I + z]/(2 + Sqrt[I + z])]], {z, -10, 10}] & ) /@ {Re, Im}; (Plot[#1[Log[(2 + I) + z - 2*Sqrt[I + z]]], {z, -10, 10}] & ) /@ {Re, Im}; (Plot[#1[Log[(2 + I) + z + 2*Sqrt[I + z]]], {z, -10, 10}] & ) /@ {Re, Im}; These plots indicate that the first function introduce the jump discontinuities in the real axis Here is the branch point of this function Solve[-2 + Sqrt[I + z] == 0, z] {{z -> 4 - I}} Before obtaing Vladimir's integral let's see a simpler example g = HoldForm[Integrate[f[z], {z - 100, 100}]] The result returned by Mathematica is wrong as we see below FullSimplify[ReleaseHold[g]] {N[%], ReleaseHold[g /. Integrate[x___] :> NIntegrate[x, MaxRecursion - > 12]]} (1/2 + I/2)*Pi + Log[((101620289 - 320288*Sqrt[10001] + 8*I*Sqrt[-21072330472601 + 1017117472601*Sqrt[10001]])/100020001)^(- (1/4) + I/4)] {1.3704665093461683 + 1.370466509346168*I, 2.941262836141067 + 2=2E941262836141067*I} However using the UNDOCUNTATED (for Integrate; it is well documentated for NIntegrate in the Help Browser. BTW, I believe that its first appearance-which I saw myself of course!-is an old post by Ted Ersek) setting of the third element in the iterator of Integrate as follows forces Mathematica to consider the above mentioned discontinuity, giving in this way the correct result! FullSimplify[Integrate[f[z], {z, -100, 4 - I, 100}]] N[%] (1/4 + I/4)*(4*Pi + Log[((101620289 - 320288*Sqrt[10001] + 8*I*Sqrt[-21072330472601 + 1017117472601*Sqrt[10001]])/100020001)^I]) 2=2E941262836141065 + 2.941262836141065*I Previous setting is equivalent to FullSimplify[Limit[F[z], z -> 100, Direction -> 1] - Limit[F[z], z -> 4 - I, Direction -> -1] + Limit[F[z], z -> 4 - I, Direction -> 1] - Limit[F[z], z -> -100, Direction -> -1]] N[%] (1/4 + I/4)*(4*Pi + Log[((101620289 - 320288*Sqrt[10001] + 8*I*Sqrt[-21072330472601 + 1017117472601*Sqrt[10001]])/100020001)^I]) 2=2E941262836141065 + 2.941262836141065*I Since Mathematica return unevaluated both Integrate[F[z], {z, - Infinity, Infinity}] and Integrate[F[z], {z, -Infinity, 4 - I, Infinity}] we will use the previous setting with Limits for our needs. We have Limit[F[z], z -> Infinity] (*correct*) 0 Limit[F[z], z -> -Infinity] (*incorrect; it misses a minus sign*) (1/2 + I/2)*Pi (Limit[F[z], z -> 4 - I, Direction -> #1] & ) /@ {-1, 1} (*correct*) {(-(1/4) - I/4)*(Pi - 2*ArcTan[1/2] - I*Log[5]), (1/4 + I/4)*(Pi + 2*ArcTan[1/2] + I*Log[5])} So, (at last!) we finally have (0 - (-(1/4) - I/4)*(Pi - 2*ArcTan[1/2] - I*Log[5])) + ((1/4 + I/ 4)*(Pi + 2*ArcTan[1/2] + I*Log[5]) - (-(1/2 + I/2))*Pi) Simplify[%] {(N[#1, 22] & )[%], NIntegrate[f[x], {x, -Infinity, Infinity}, WorkingPrecision -> 50, PrecisionGoal -> 20]} (1/2 + I/2)*Pi + (1/4 + I/4)*(Pi - 2*ArcTan[1/2] - I*Log[5]) + (1/4 + I/4)*(Pi + 2*ArcTan[1/2] + I*Log[5]) (1 + I)*Pi {3.1415926535897932384626433832795028842`22. + 3=2E1415926535897932384626433832795028842`22.*I, 3.14159265358979323846264338327950287583`21.64106815259695 + 3=2E14159265358979323846264303665222933418`21.64106815259695*I} I=2Ee the desired integral is equal to (1 + I)*Pi. Best Regards Dimitris =CF/=C7 Michael Weyrauch =DD=E3=F1=E1=F8=E5: > Hello, > > Dimitris, this is a marvelous solution to my problem. I really apprecia= te > your help. I will now see if I can solve all my other (similar) integrals= using the same trick. > Timing is not really the big issue if I get results in a reasonable amoun= t of time. > > Also the references you cited are quite interesting, because they give so= me insight > what might go on inside Mathematica concerning integration. > > I also agree with you that it is my task as a programmer to "help" Mathem= atica > and guide it into the direction I want it to go. Sometimes this is an "ar= t"... > > Nevertheless, in the present situation I do not really understand why Mat= hematica wants > me to do that rather trivial variable transformation, which is at the hea= rt of your solution. > The integrand is still a rather complicated rational function of the same= order. The form > of the integrand did not really change > substantially as it is the case with some other ingenious substitutions o= ne uses in order to > do some complicated looking integrals "by hand". > > I think the fact that we are forced to such tricks shows that the Mathema= tica integrator > is still a bit "immature" in special cases, as also the very interesting = article by D. Lichtblau, > which you cite, seems to indicate. So all this is probably perfectly know= n to the > Mathematica devellopers. And I hope the next Mathematica version has all = this "ironed out"?? > > Many thanks again, Michael