Monte Carlo calculation of volumes of n-balls

[dcromley]

I could have sworn that somebody started a thread here regarding an article in The American Scientist.

http://www.americanscientist.org/issues/pub/an-adventure-in-the-nth-dimension/1

So I looked at it, got interested, and wanted to verify some of the "zaniness" of the subject. He is working with "n-balls" of radius 1, which are "spheres" in n dimensions. A 3-ball is the normal sphere in 3 dimensions; the volume is 4/3pi*r^2 = 4.1888. A "2-ball" is a circle in 2 dimensions; the "volume" is 2*pi*r = 3.1416. A 1-ball is a line from -1 to 1; the "volume" is 2. Going the other way -- going up in dimensions is not clear at all, but he gives a table on page 4 and a formula on page 8.

My thought was to write a program to divide the "n-cube" into 0.01 increments and count the points that were "inside" the "n-ball" by comparing the distance from the center to 1. Oops -- in 10 dimensions, that would be 100^10 = 100,000,000,000,000,000,000 points. Maybe increments of 0.1. That's still 10^10 = 10,000,000,000 points.

Another approach is the "Monte Carlo" method

http://en.wikipedia.org/wiki/Monte_Carlo_method

which is picking many random points and counting the ones that are inside the "sphere". The distance from the center is the square root of the sum of the squares of the n coordinates. The volume is the ratio of points "inside" to total points times the volume of the "n-cube".

This TB program does that, and surprise, surprise -- the results agree with the table! The higher dimension "n-balls" are slow in converging because the "volume" of the ball is very small compared to the vast "volume" of the "n-cube" of size 2. For dimension 10, the "n-ball" is 2.55 and the "n-cube" is 1024 cubic units.

Changing the program to size 1 (n-ball in a "unit" n-cube") gives more normal results -- the volumes don't increase.

Table of volumes of n-balls
Dim  Size 2  Size 1
1    2.0000  1.0000 
2    3.1416  0.7854 
3    4.1888  0.5236 
4    4.9348  0.3084 
5    5.2638  0.1645 
6    5.1677  0.0807 
7    4.7248  0.0369 
8    4.0587  0.0159 
9    3.2985  0.0064 
10   2.5502  0.0025

  Uses "Console", "UI"
  Global hWnd As Long
  hWnd = Canvas_Window "    " & _ 
    "Monte Carlo Calculation of Volumes of ""N-Balls"" in size 2 ""N-Cubes"" in Dimensions 1-10;" & _ 
    "    ESC To Exit", 0, 0, 1018, 736
  ' x goes from 0-1017; y from 0-735
  ' dCromley, Dec, 2011
  Const ymax As Number = 6 ' for y scaling
  Const maxdim As Number = 10 ' max dimensions
  Canvas_Attach hWnd, 0, %TRUE
  Canvas_Scale(0,ymax, 1017,0)
  Canvas_Clear Rgb(255,255,255)
  Global aIns(maxdim), aIters(maxdim) As Long, aVol(maxdim) As Single
  Global canvasx, niterations As Long ' ix is screenx

Sub TBMain()
' http://www.americanscientist.org/issues/pub/an-adventure-in-the-nth-dimension/
  Local i1, i2, idim As Long, rsqr, vol As Single
  Randomize
  ygrid()
  For canvasx = 1 To 1017
    niterations = canvasx^2 ' 1 to 1000000 or so
    For idim = 1 To maxdim ' for each of 10 dimensions
      vol = GetVol(idim) ' volume is ratio of ins / ins+outs
      vol *= 2^idim ' adjust for size 2 sides
      aVol(idim) = vol ' save for printing
      Canvas_SetPixel(canvasx, vol)      
    Next idim                   
    If canvasx > 100 And Mod(canvasx, 10) = 0 Then zPrint() ' every once in awhile
    If GetAsyncKeyState(%VK_ESCAPE) Then Exit Sub
  Next canvasx         
  Beep
  Canvas_WaitKey
End Sub                         

Function GetVol(nDim As Long) As Single
' -- get random points in a n-dimension box 
'    return vol, the fraction of points within ball
  Local Ins, iIter, iDim As Long, rsqr, vol As Single
  For iIter = aIters(ndim)+1 To nIterations
    Incr aIters(ndim) 
    rsqr = 0
    For iDim = 1 To nDim ' add up the iDim sides^2
      rsqr += (Rnd-.5)^2 ' -.5 to +.5 squared
    Next iDim           
    ' don't take sqrroot, just comare with .25
    If rSqr < .25 Then Incr aIns(ndim) ' if within ball, incr Ins
  Next iIter                             
  vol = aIns(ndim)/aIters(ndim) ' volume is fraction of Ins / Ins+Outs
  Function = vol
End Function

Sub zPrint() ' labels and values
  Local idim, ix As Long
  For idim = 1 To maxdim     
    ix = canvasx - idim*10
    If Mod(ix, 100) = 0 Then 
      Canvas_SetPos(canvasx, aVol(idim))
      Canvas_Print idim ' label
    End If  
    Canvas_SetPos(80+idim*80, .8)
    Canvas_Print idim & ": " & Format$(aVol(idim), "0.0000")
  Next idim                     
  Canvas_SetPos(50, .8)
  Canvas_Print "Iters="&niterations ' # iterations
  Canvas_Redraw
End Sub

Sub ygrid()
  Local iy As Long
  For iy = 0 To ymax-1
    Canvas_SetPos(0, iy+.1)
    Canvas_Print Format$(iy, "0.00")
    Canvas_Line((0, iy), (1017, iy), Rgb(192,192,192))
  Next iy
End Sub

[danbaron]

Here is a similar problem which I solved with thinBasic, using the Monte Carlo method.

Dan

http://www.thinbasic.com/community/s...8473#post78473

[dcromley]

Very interesting also, thanks.

As I was running mine, I knew how dependent it was on "rnd". And reading your and Eros' dialog on random number generation was interesting. Quite a bit of analysis to that. I'll look for "Numerical Recipes in C". Regards, Dave