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