Great Michael Trott & six-fold integral

*To*: mathgroup at smc.vnet.net*Subject*: [mg30124] Great Michael Trott & six-fold integral*From*: seidovzf at yahoo.com (Zakir F. Seidov)*Date*: Sun, 29 Jul 2001 21:26:16 -0400 (EDT)*Organization*: The Math Forum*References*: <9jtk3e$87k$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

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