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