rss

Mathematical digression : Buffon's needle

5

Category : C#, English, General, Programming

Buffon’s needle is a statistic experiment created by Comte George-Louis Leclerc (sounds so frenchie).

The principle is to drop needles on a parquet floor and check if the needle cross one of the parquet line (providing each parquet’s strip has the same height).

The experiment can be modeled as two random variables representing, for the first, the distance between the center of the needle and the closest parquet line, for the second, the acute angle between the needle and the line. Both random variables follow a uniform distribution, respectively of parameters (0, {\frac{length(needle)}{2}} ) and parameters (0, {\frac{\pi}{2}} ).

A needle cross a line when the angle is superior to 0 and, depending on the previous angle, the distance is less than :

{sin(angle)\times\frac{height(strip)}{2}}

Thus, the condition « the needle cross the line » can be summarized by this equation :

{x \leq \frac{l}{2} \sin (\theta)}

Where a is the length of the needle, θ the angle and x the distance.

The cool part is that computing the joint probability of the two random variable with the crossing condition gives you an expression of π :

{\pi = \frac{2aN}{ln}}

Where N is the number of needle dropped, l the height of a strip and n the number of needle that crossed the line.

Therefore, for a good number of needle drop you can get a value fairly close to π.

The following C# source code is doing precisely that :

using System;

namespace BuffonApp
{
  class BuffonNeedles
  {
	public static void Main ()
	{
	  // Number of time we will drop a needle
	  const int N = 1000000;
	  // Interval between each line
	  const int a = 10;
	  // Size of a needle
	  const int l = a - 4;
	  // Our random number generator
	  Random r = new Random ();
	  // Number of time we cross the line
	  int n = 0;

	  for (int i = 0; i < N; i++) {
		double x = r.NextDouble () * (a / 2);
		double theta = r.NextDouble () * (Math.PI / 2);
		if (x <= (l / 2) * Math.Sin (theta))
		  n++;
	  }

	  Console.WriteLine ("Pi approximation : " + ((double)(N * 2 * l)) / (n * a));
	}
  }
}

Which give us a rather good approximate of : {\pi \approx 3,14208745948547}

Comments (5)

This method is also referred to as the Monte Carlo method, although it’s just one application among others in the Monte Carlo algorithm class.

using the sine function to approximate Pi, ss a chicken egg problem for me :)

Of course, if you have access to a Sin function, you already have an easier way to get to Pi – just find out what x makes Sin(x) = 1, and double it. And of course, your sample code is using PI directly, albeit to range-limit an angle selection.

A more interesting method, in a didactic sense, is to estimate the proportion of the area of a circle to the area of its smallest enclosing square. The area of the circle is pi*r*r, while the area of the square is 4*r*r. Dividing one by the other, as estimated by needle drops (x,y) where both x and y are in the domain -r..r and tested with circle equation of x*x+y*y=r*r (>r*r => outside circle), and you have a way of probabilistically estimating pi without any hidden pi in the methodology.

Here’s some sample code:

using System;

class App
{
const int Iter = 1000000;

static void Main()
{
Random r = new Random();
double radius = 0.5;
double radiusSquared = radius * radius;
int hitCount = 0;

for (int i = 0; i < Iter; ++i)
{
double x = r.NextDouble() * radius * 2 – radius;
double y = r.NextDouble() * radius * 2 – radius;
if (x * x + y * y <= radiusSquared)
++hitCount;
}

// area of circle / area of square
// = pi * r * r / 4 * r * r
// = pi / 4
double proportion = hitCount / (double) Iter;
Console.WriteLine(« Pi estimated at {0} », proportion * 4);
}
}

@gliss:
It’s yet another one of those « repeat it a lot to see where it converges » :-) .

@shamaz:
Yeah, I coded it just for fun (and because I don’t have the luxury to throw needles on my parquet all day long :-) ).

@Barry:
Nice one.

I don’t think there is really a « hidden » pi in the method since you can express your sin function without any kind of pi (via series for instance). But, agreed, it’s disputable :-) .

Amusant, il y a quelques semaines, au cours de math, nous avons vu les intégrales (ouai bon je suis toujours en secondaire :P ) et j’ai trouvé un moyen très rapide d’avoir une approxymation de PI (sans avoir un accès aux fonctions trigonometriques). Testé sur une Casio en CasioBasic \o/.

Les intégralles permettent d’obtenir l’aire d’une fonction. Et PI est le rapport entre la circonférence du cercle et son rayon, mais nous savons aussi que PI*R² = Aire.
Il suffit donc de prendre la formule du cercle (en réalité le demi-cercle pusique une fonction ne connais que une image pour un X): Y = SQRT(R-X²).
Et d’y calculer son aire entre 0 et 1 pour un demi-cercle de 1 de rayon et de la multiplier par 4 (puisque 0 et 1 ne sont que un quart de cercle), de là nous aurons directement PI (Si R = 1 => PI*R² = PI, PI est donc l’aire du cercle). Les intégrales nous permettent de calculer cette aire.

On peu avoir une estimation d’une intégralle en additionant des échantillon d’aire de la fonction en « testant » la fonction sur plusieurs valeurs de X. Plus nous testons d’échantillons, plus la valeur sera précise : http://files.getwebb.org/index.php?mode=view&id=cre66i0p

Et le code: http://pastebin.com/f447dec0c

On peut voir que l’estimation est beaucoup plus proche de PI, on peut atteindre le même degré de précision avec seulement 100 itérations :-) . Mais ca reste assez lourd à cause de SQRT() (même si 100 itération doivent quand même s’exécuter plus vite que 1000000).

Il existe aussi une fonction beaucoup plus simple mais il ne s’agit pas d’une approximation et le temps de calcul de chaque décimale et relativement exponentielle:
PI = 1 / 1 – 1 / 2 + 1 / 3 …((-1 * (n+1)%2) * 1) / n

Post a comment