Date: Sat, 27 Jan 96 17:18:27 CST
From: rusin (Dave Rusin)
To: borrisd@magi.com
Subject: Re: spheric
You asked,
>What is the average distance between 2 randomly selected points
>INSIDE a sphere (3D) of radius R (uniform distribution)?
>
>I'm fairly confident this is the way to solve this; but being presently
>quite rusty in this area, I can't work it out.
>I'd appreciate the answer and how you worked it out. Thank you.
>
>[BALL (sphere) of radius R = B(R)].
>Pick an arbitrary random point (a,b,c) in B(R). Its distance to any
>other point (x,y,z) in B(R) is the square root of the function:
>f(a,b,c,x,y,z) = (x-a)^2 + (y-b)^2 + (z-c)^2.
>Next, to get the average distance between (a,b,c) and points randomly
>picked in B(R), integrate the sq. root of this function over B(R) against
>dxdydz and then divide by the volume of B(R). (The volume is 4/3 pi R^3.)
>This gives a new function F(a,b,c). Since (a,b,c) was randomly picked,
>integrate F over B(R) and divide again by the volume.
Exactly so.
The only problem seems to be doing the integrals. Of course we expect
an answer of the form c*R where c is an absolute constant, greater
than zero, and no more than 2. My calculations give (453 / 280 ) R.
I guess I would have expected something a little smaller, but not much,
so I give this answer with some confidence.
To do the integrals, proceed like this. For a fixed (a,b,c), the
average of the distance function will be the same as the average which
results from any other (a', b', c') lying on the same sphere
centered at the origin, because of the rotational symmetry of the ball.
So to compute this average it will suffice to do the computation
using the point (a', b', c') = (R1, 0, 0), where R1^2=a^2+b^2+c^2.
(so 0 <= R1 <= R)
All right then, now we have to find the average distance from (R1,0,0).
As you suggest we should integrate over the whole ball of (x,y,z)'s.
Do it in stages; first integrate over each vertical disc x=const. You can
do that integral in polar coordinates, since the points on that disc
are (x,y,z)=(const, r cos v, r sin v) with 0 <= v <= 2 pi and
0 <= r <= sqrt(R^2 - x^2). So we need the integral
sqrt(R^2-x^2) 2 pi
Int Int sqrt(r^2 + (R1-x)^2) r dv dr
0 0
The integrand doesn't depend on v so that integral just contributes
a factor of 2 pi. The other integral we can do substituting
u = ( r^2 + (R1-x)^2 ) so that du = 2 r dr and the integral is
u=(R^2-x^2)+(R1-x)^2
2 pi Int sqrt( u ) (1/2)du
u=(R1-x)^2
The indefinite integral is (1/3)(u^(3/2)). Plugging in the endpoints gives
the value
2 pi (1/3) [ R^2 + R1^2 - 2 R1 x ] ^ (3/2) - 2 pi (1/3) |R1-x|^3
This is the integral of the distance function over the disc.
This is what we're supposed to integrate from x= -R to R.
Do this by another substitution if you like, but it's a linear function
in x to the 3/2 power, minus a polynomial in x times a sign function,
so we quickly get
2 pi (1/3) * (2/5) [ R^2 + R1^2 - 2 R1 x ] ^(5/2) / [-2 R1]
-2 pi (1/3) * (-1/4) (R1-x)^4 * sign(R1-x)
This time plugging in the endpoints gives
(2/15) pi / R1 * [ |R + R1|^5 - |R - R1|^5 ]
-(1/6) pi * [ (R1-R)^4 + (R1+R)^4 ]
If you divide this by the volume of the ball (= (4/3) pi R^3 ) you get
the average distance from the point (R1, 0, 0) to all the other points
in the ball of radius R.
Now you want to let the other point (a,b,c) (considered fixed in the
previous discussion) vary over the ball. It's best to do this one
with spherical coordinates, since as we have noted the integral of the
distances from (a,b,c) to the other (x,y,z) depends only on the
distance R1 from (a,b,c) to the origin. So when we integrate over
a fixed sphere, we get the previous integral, multiplied by the surface
area 4 pi R1^2 of that sphere:
4 pi R1^2 * { (2/15) pi / R1 * [ (R + R1)^5 - (R - R1)^5 ]
-(1/6) pi * [ (R1-R)^4 + (R1+R)^4 ] }
= (16/15) pi^2 R1 * [ 5 R1 R^4 + 15 R1^3 R^2 + R1^5 ]
- (4/3) pi^2 R1^2 * [ R^4 + 6 R1^2 R^2 + R1^4 ]
= pi^2 R1^2 * [ 4 R^4 + 8 R1^2 R^2 - (4/15) R1^4 ],
if I did my algebra right.
This then is integrated over the interval [0, R]. Again, for polynomials,
no problemo:
pi^2 * R^7 * [ 4/3 + 8/5 - 4/105 ] = 304/105 pi^2 R^7.
OK, so divide this twice by (4/3) pi R^3 to get the average distance.
I get (453 / 280 ) R.
You might want to double check my arithmetic; perhaps if I get a chance I'll
have Maple check the integrals.
A more interesting problem than I would have thought!
dave
==============================================================================
Date: Sun, 28 Jan 96 00:49:21 CST
From: rusin (Dave Rusin)
To: borrisd@magi.com
Subject: Re: spheric
Something bugged me about my answer; turned out it was off a little.
Recall that you had asked,
>What is the average distance between 2 randomly selected points
>INSIDE a sphere (3D) of radius R (uniform distribution)?
I computed (453/280)R; this is wrong. I began with this problem:
>All right then, now we have to find the average distance from (R1,0,0).
>As you suggest we should integrate over the whole ball of (x,y,z)'s.
>Do it in stages; first integrate over each vertical disc x=const.
This first stage I wrote as
> sqrt(R^2-x^2) 2 pi
> Int Int sqrt(r^2 + (R1-x)^2) r dv dr
> 0 0
It in turn had to be integrated over x from -R to R. I've now done all the
calculations in Maple, and gotten the same answer as I got by hand, although
Maple chose to simplify it at this stage and so expressed the answer as:
2 2 4 4
1/15 Pi (10 R R1 - R1 + 15 R )
Next I explained why I had to integrate this times 4 pi R1^2, an expression
which I simplified to
> = pi^2 R1^2 * [ 4 R^4 + 8 R1^2 R^2 - (4/15) R1^4 ],
but as can be seen from Maple's output,this is a goof; the "8" should
be 8/3. Of course this changes the integral of this over [0,1]; Maple
says this is 64/35 * pi^2 * R^7 (I had the fraction as 304/105)
>OK, so divide this twice by (4/3) pi R^3 to get the average distance.
Now this is (36/35) R. (Maple agrees.)
Uncertain as to whether to trust my algebra (now proven once wrong) or
my intuition (now corroborated -- I thought the other number was too big!)
I ran a simulation instead. Here's a UBASIC program:
10 randomize
20 dim X(6)
30 S=0:N=0
40 for I=1 to 100000
50 for J=1 to 6
60 X(J)=2*rnd-1
70 next J
80 if X(1)^2+X(2)^2+X(3)^2>1 then 130
90 if X(4)^2+X(5)^2+X(6)^2>1 then 130
100 D=sqrt((X(1)-X(4))^2+(X(2)-X(5))^2+(X(3)-X(6))^2)
110 S=S+D
120 N=N+1
125 'next four lines for pretty-printing the answer
130 if I@1000>0 then 160
140 print I,N,S/N,(S/N)*35
145 if I@10000=0 then 160
150 locate ,posy-1
160 next I
I ran a couple of trials of 10000 pairs of points; one showed the mean
distance to be just over 36/35; the other, just under. I'm satisfied.
dave