MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

Great Michael Trott & six-fold integral


Yes, Michael,

Your cunning usage of Mathematica is really amazing.

I'm not sure that ordinary user (as me, e.g.) 
could understand your session, but I'd try...

I've just rerun your session and  got result in seconds!!

Unfortunately (but as usually) I've made a lot of misprints:
missing closing ")", sign and factor 2...

Thank you very much indeed for your great posting,

But I dare ask:
if you can, with your method, also solve 
the general case of  arbitrary limits 

Int_0 ^a Int_0 ^b Int_0 ^c Int_0^ A Int_0 ^B Int_0 ^C  

which was solved in my preprint, please?!


Zakir F. Seidov

seidov at bgumail.bgu.ac.il

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Subject: [mg30124]      Re: six-fold integral
Author:       Michael Trott <mtrott at wolfram.com>
Organization: Steven M. Christensen and Associates, Inc and
MathTensor, Inc.

Your result 

(2 3^(1/2)-2^(1/2)-1)/5 + pi/3 + ln[(2^(1/2)-1)(2-3^(1/2)]

misses a closing ')'. Assuming you mean 

(2 3^(1/2)-2^(1/2)-1)/5+Pi/3 + Log[(2^(1/2)-1)(2-3^(1/2))] // N,

this is unfortunately negative (-0.941156).

Using the joint probability density for (x - y) it is straightforward
to calculate the following exact result in just a few lines:

In[1]:= logRewrite[expr_] := FixedPoint[PowerExpand[# //.
         Log[a_] :> Log[Factor[Together[a]]]]&, TrigToExp[expr]]

In[2]:= indefiniteIntegral = Integrate[# 2(1 - z), z]& /@
Expand[logRewrite[
   Integrate[Integrate[1/Sqrt[x^2 + y^2 + z^2] 2(1 - x), x] 2(1 - y),
y]]];

In[3]:= replaceLimits[expr_, xyz_] :=((# /. xyz -> 1) - (# /. xyz ->
0))&[
 Expand[logRewrite[expr]] /. xyz^n_?Positive Log[xyz] -> 0]

In[4]:= Fold[replaceLimits, indefiniteIntegral, {z, y, x}] //
Together;

In[5]:= FullSimplify[%]

Out[5]= -((2*Pi)/3) + Log[3 + 2*Sqrt[2]] + 
  1/10*(4 + 4*Sqrt[2] - 8*Sqrt[3] - 25*Log[2] + 50*Log[1 + Sqrt[3]] - 
    5*Log[2 + Sqrt[3]])

In[6]:= N[%]

Out[6]= 1.88231

This agrees with a Monte-Carlo check:

In[7]:=
SeedRandom[888];
Compile[{{n, _Integer}},
        Module[{sum = 0.},
             Do[sum = sum + 1/Sqrt[(Random[] - Random[])^2 +
                                   (Random[] - Random[])^2 +
                                   (Random[] - Random[])^2],
                {n}]; sum/n]][10^6]

Out[7]= 1.88295




  • Prev by Date: Re: Graphics from Mathematica to MS Word
  • Next by Date: Two factors of (10^71-1)/9 = R71
  • Previous by thread: Re: Graphics from Mathematica to MS Word
  • Next by thread: Two factors of (10^71-1)/9 = R71