# Thread: Monte Carlo calculation of volumes of n-balls

1. ## Monte Carlo calculation of volumes of n-balls

I could have sworn that somebody started a thread here regarding an article in The American Scientist.
HTML Code:
`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
HTML Code:
`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()
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
```

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

Dan

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

3. 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

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•