1
$\begingroup$

I am trying to numerically reproduce the Mathematica result for the density of strings in a sphere as explained here.

  1. The program is to generate random, uniformly distributed pair points on a unit sphere. This pair forms a line segment.

  2. Use shells of constant thickness, incrementally increasing in radius, to calculate the density of lines at that radius by finding the portion of all line segments that lie in the shell.

  3. Plot the density as a function of r.

Code:

noOfTrials = 500;
points = RandomPoint[Sphere[], {noOfTrials, 2}];
(*gives the length of a line from pt1 to pt2 that lies inside a 
 sphere at origin of radius r*)
dFunc[pt1_, pt2_, r_] := Module[
  {midPt, h},  
  midPt = (pt1 + pt2)*0.5;
  h = Norm[midPt];
  If[h >= r, 0, 2 N@Sqrt[r^2 - h^2]]
]  
(*simulation*)
sim = Module[
{dr, l, dlTotal},
dr = 10^-2;
 Table[
 dlTotal = 
  Plus @@ Table[(dFunc[#1, #2, r + dr] - dFunc[#1, #2, r]) & @@ 
     pair, {pair, points}];
 {r, dlTotal/(4 Pi r^2 dr)}
 , {r, 0.01, 1, dr}
 ]
];  
(*plot*)
ListPlot[sim]

Problem

  1. Even for up to 10000 pts, the plot doesn't seem to tend to a constant value near 1/$\pi$ as it should as seen from the answer mentioned before

I've separately checked and rechecked the dFunc, which I presumed to be the only source of error, but found it to be correct; tweaked dr value; Checked the geometry -- but nothing seems wrong.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
0
$\begingroup$

Update after comments: It turned out that the problem with OP's code was not dividing by the number of lines.

Here is a straightforward, alternative way to compute the same thing:

distance[Point[{pt1_, pt2_}]] := Norm[pt1 - pt2]
distance[_] := 0

massInsideSphere[line_, r_] := 
 distance@RegionIntersection[Sphere[{0, 0, 0}, r], line]

density[lines_, r_] := Module[{mass1, volume1, mass2, volume2},
  mass1 = Total[massInsideSphere[#, r] & /@ lines];
  volume1 = Volume[Ball[{0, 0, 0}, r]];
  mass2 = Total[massInsideSphere[#, r + 0.01] & /@ lines];
  volume2 = Volume[Ball[{0, 0, 0}, r + 0.01]];
  (mass2 - mass1)/((volume2 - volume1) Length[lines])
  ]

nrOfTrials = 10000;
lines = InfiniteLine /@ RandomPoint[Sphere[], {nrOfTrials, 2}];

estimates = Table[{r, density[lines, r]}, {r, 0.1, 0.9, 0.1}];
ListPlot[estimates]

Mathematica graphics

This agrees approximately with the theoretical result $\rho = \frac{1}{\pi}$.

1./Pi

0.31831

Improve this answer
$\endgroup$
3
  • $\begingroup$ 1. The perpendicular bisector of a chord always passes through the center. 2. I have used (dFunc[#1, #2, r + dr] - dFunc[#1, #2, r]) - is that not the length of line in the shell? 3. Your density function computes the total mass in sphere/sphere volume. How is that the local density at r? Shouldn't it be Total mass in shell/shell volume? I think you have computed the average density of spheres of different r instead of the radial variation of the local density in a unit sphere? $\endgroup$
    – lineage
    Apr 28, 2019 at 17:39
  • $\begingroup$ @lineage ok, you are right about your answer. About my answer, sure, if you want to make sure that the density locally around some $r$ is $1/\pi$ then you have to do the mass computation twice for two different $r$ close to each other, and then divide the difference by the volume between those radii. It's easy to do by modifying the code that I've given. $\endgroup$
    – C. E.
    Apr 28, 2019 at 18:07
  • $\begingroup$ thanks to you i realised that local density will be proportional to the number of lines,turns out it was indeed the only mistake(note that division by dr is still necessary)--please remove the misleading code part of your answer and i shall accept it--thanks again for the help. $\endgroup$
    – lineage
    Apr 28, 2019 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.