Re: six-fold integral
- To: mathgroup at smc.vnet.net
- Subject: [mg30098] Re: six-fold integral
- From: mtrott at wolfram.com (Michael Trott)
- Date: Sat, 28 Jul 2001 01:51:00 -0400 (EDT)
- References: <9jr6uk$4mn$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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