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